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1.  INTRODUCTION 

Under  this  contract,  Northwest  Research  Associates  (NWRA)  is  performing  measurements 
of  the  ionosphere  in  order  to  investigate  plasma  phenomena  in  both  the  natural  ionosphere  and  in 
the  ionosphere  perturbed  by  the  High-frequency  Active  Auroral  Research  Program  (HAARP) 
high-frequency  heater.  The  research  being  conducted  falls  within  Hanscom  Technical  Areas 
3(a),  “Ionospheric  Effects  Research  and  Department  of  Defense  Systems,”  and  3(b), 
“Ionospheric  Research  Technology,”  of  the  Broad  Agency  Announcement  VS-03-01  released  by 
the  Air  Force  Research  Laboratory  (AFRL). 

As  a  contribution  to  Technical  Area  3(b),  we  collaborated  with  scientists  and  engineers  from 
AFRL,  the  Naval  Research  Laboratory  (NRL),  and  other  research  organizations  in  the 
application  of  diagnostic  instrumentation  to  HAARP.  Among  the  HAARP  instruments  that  we 
applied  to  ionospheric  research  are  four  NWRA  ITS  1  OS  coherent  radio  receiving  systems  for 
measuring  relative  Total  Electron  Content  (TEC)  and  recording  scintillation,  a  mobile  NWRA 
ITS30m  coherent  radio  receiving  system  for  recording  TEC  and  scintillation  in  support  of 
HAARP  campaign  activities,  and  an  Ashtech  Model  Z-FX  GPS  receiver  for  measuring  absolute 
TEC.  TEC  from  these  instruments,  and  phase-scintillation  records  from  the  ITS  1  OS  receivers, 
are  posted  on  the  HAARP  web  site  for  telescience  applications  and  for  decision-making  during 
active  experiments,  and  the  TEC  data  were  inverted  using  NWRA’s  tomography  software  to 
produce  images  of  the  F  layer  over  Alaska,  with  the  resulting  images  posted  on  the  HAARP  web 
site.  Research  topics  under  the  foregoing  objectives  are  reported  in  Section  2  and  Appendixes  A, 
B,  and  C. 

In  Technical  Area  3(a),  topics  of  research  were  (1)  studies  pertinent  to  improvement  of 
ionospheric  scintillation  models,  (2)  various  tasks  focused  on  developing  techniques  for 
partitioning  ionospheric  TEC  derived  from  GPS  measurements  into  ionospheric  and 
plasmaspheric  content  and  employing  them  to  study  the  behavior  of  the  plasmapause  and 
plasmaspheric  TEC  using  GPS  measurements  of  TEC,  and  (3)  studies  of  the  sources  and 
magnitude  of  uncertainty  in  TEC  derived  from  GPS  measurements.  Research  topics  under  the 
foregoing  objectives  are  reported  in  Section  3. 

Note  that  work  accomplished  during  the  period  17  February  2004  through  17  August  2005 
were  covered  by  Scientific  Report  No.  1  ( Secan  et  al.,  2005).  Materials  presented  in  that  report 
will  not  be  repeated  here  with  the  exception  of  changes  that  have  been  made  during  the  work 
completed  since  the  end  of  that  report.  We  will  refer  to  that  report  as  SRI  in  this  document. 

2.  HAARP  TOPICS  -  TECHNICAL  AREA  3(B) 

Under  this  contract,  NWRA  coordinated  installation  of  enhanced  diagnostic  instrumentation 
and  participated  in  research  employing  some  of  them.  Our  activities  aimed  at  enhancing 
HAARP  diagnostics  were  carried  out  primarily  by  NWRA  consultants  A.  Lee  Snyder,  Edward  J. 
Kennedy,  Jens  Ostergaard,  and  Spencer  Kuo.  NWRA  also  employed  two  summer  student 
interns,  Mr.  Troy  Lawlor  and  Ms.  Hayley  Crow,  who  worked  with  AFRL  personnel  at  Hanscom 
AFB  and  at  the  HAARP  site. 
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2.1.  Summary  of  HA  ARP  Data  Collection  Activities 

NWRA‘s  HAARP  research  activities  centered  on  collection  and  analysis  of  transionospheric 
radiowave  data.  Transmission  of  two  or  more  phase-coherent  radio  signals  from  satellites  above 
the  ionosphere  to  receivers  on  the  ground  permits  measurement  of  the  path  integral  of  plasma 
density  (so-called  total  electron  content,  or  TEC)  through  the  ionosphere.  TEC  is  proportional  to 
the  dispersive  (differential)  phase  between  the  two  signals  (and,  with  less  precision  but  greater 
certainty,  to  differential  group  delay).  Recording  of  rapid  fluctuations  in  dispersive  phase  and  in 
received  signal  strength  yields  measurement  of  trans-ionosphere  radiowave  scintillation. 

Under  this  contract,  NWRA  measured  TEC  on  slowly  moving  paths  penetrating  the  entire 
ionosphere  by  recording  dispersive  phase  and  differential  group  delay  registered  on  signals 
transmitted  from  satellites  of  the  Global  Positioning  System  (GPS)  to  Gakona,  AK.  We  also 
measured  both  TEC  and  scintillation  on  paths  between  satellites  moving  rapidly  over  Alaska  in 
high-inclination,  low-earth  orbits  (LEO)  just  above  the  main  ionospheric  (F)  layer  and  several 
ground  stations.  Collection  and  analysis  of  the  LEO  data  sets  is  presented  in  Section  2.1.1.,  and 
collection  and  analysis  of  the  GPS  data  sets  is  presented  in  Section  2. 1 .2.  All  of  the  data  and 
displays  described  in  these  two  sections  are  archived  on  either  the  HAARP  computer  network  at 
the  HAARP  facility  near  Gakona,  AK,  or  on  servers  at  NWRA’s  Bellevue  office. 

2.1.1.  LEO  Satellite  and  Tomography  Data  Sets 

NWRA  operates  NWRA  ITS  1  OS  receiver  systems  at  four  sites  in  Alaska  in  support  of 
HAARP  operations:  Cordova,  Gakona  (at  the  HAARP  facility),  Delta  Junction,  and  Arctic 
Village.  In  addition,  we  obtain  data  from  two  University  of  Texas  at  Austin  (UTA)  CIDR 
receivers  operated  by  the  University  of  Alaska  Geophysical  Institute  (UAF  GI)  at  two  sites  in 
Alaska:  Ft.  Yukon  and  Kaktovik.  Figure  1  shows  the  locations  of  these  sites  on  a  map  of 
Alaska.  The  site  at  Arctic  Village  was  established  during  the  current  contract,  and  is  an 
important  additional  site  in  the  tomography  chain.  Data  from  LEO  satellite  passes  collected  at 
NWRA’s  four  sites  are  all  displayed  in  this  format  on  the  HAARP  web  site,  accessible  through 
that  site’s  Data  Index  page. 


Figure  1.  Location  of  LEO  receiver  sites  in  Alaska.  Sites  indicated  by  stars  (C:  Cordova,  G:  Gakona,  D: 
Delta,  A:  Artie  Village)  have  NWRA  ITS10S  receivers.  Sites  indicated  by  circles  (Y:  Ft.  Yukon,  K: 
Kaktovic)  have  UTA  CIDR  receivers. 
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Data  from  these  sites  are  transferred  to  both  HAARP  facility  computers  and  to  NWRA’s 
Bellevue  office  for  generation  of  displays  as  shown  in  Figure  A1  and  tomographic  analyses. 
Descriptions  of  the  tomographic  processor  used  to  generate  images  and  the  criteria  used  to  select 
passes  for  tomographic  analysis  can  be  found  in  Andreasen  et  al.  (2004)  and  Fremouw  et  al. 
(1992,  1994).  The  implementation  of  the  processing  algorithms  and  the  associated  data 
preprocessing  was  described  in  detail  in  SRI,  and  examples  of  the  various  products  (data  files 
and  displays)  generated  during  the  processing  were  also  described  in  that  document. 

Under  the  present  contract,  we  produced  7,772  images  from  LEO  passes  over  the  HAARP 
tomography  chain  covering  the  interval  15  February  2004  through  31  July  2008,  all  of  which 
were  posted  on  the  HAARP  web  site.  Tables  summarizing  the  image  coverage  over  this  period 
can  be  found  in  Appendix  A. 

One  change  of  significance  that  has  been  made  to  the  tomography  analyses  posted  on  the 
HAARP  Web  site  is  the  inclusion  of  (a)  estimates  of  uncertainty  in  the  plots  of  TEC  and  f0F2 
derived  from  the  tomography  images,  and  (b)  observations  of  f0F2  from  the  Gakona  and  College 
Digisondes  when  they  are  available  at  times  near  the  image  time.  The  uncertainty  bars  are  based 
on  the  algorithm  described  later  in  this  report.  Figure  2  is  an  example  of  the  tomography  results 
as  they  are  now  plotted  on  the  HAARP  web  site  showing  these  new  additions.  The  asterisks  on 
the  f0F2  plot  section  of  the  figure  are  from  the  Gakona  and  College  Digisonde  (from  left  to  right 
on  the  figure).  The  times  given  are  the  times  the  ionograms  were  made  for  the  different  f0F2 
estimates. 

2.1.2.  GPS  Data  Sets 

NWRA  operated  an  Ashtech  Z-FX  Continuously  Operating  Reference  Station  (CORS), 
consisting  of  a  12-channel  GPS  receiver  and  an  antenna,  at  the  HAARP  Gakona  facility  for  the 
entire  period  covered  by  this  report.  Data  collection  is  performed  by  means  of  the  NWRA  GPS 
Ionospheric  Observation  System  (GIOS)  program,  supplemented  by  a  real-time  process  to 
convert  raw  data  records  from  the  GIOS  software  into  calibrated  estimates  of  equivalent  vertical 
TEC  (VTEC)  using  the  SCORE  (Self-Calibration  of  Range  Errors)  process  developed  by  NWRA 
for  the  USAF  Ionospheric  Measuring  System  (IMS)  instruments  {Bishop  et  al.,  1995).  These 
data  are  stored  in  a  database  and  are  used  to  generate  plots  of  absolute  TEC  displayed  on  the 
HAARP  Web  site.  Examples  of  the  products  generated  in  this  processing  were  given  in  SRI. 

In  addition  to  the  Ashtech  Z-FX  system,  NWRA  also  operates  a  GPS-based  CNS  clock  as  a 
time  standard  for  the  NWRA  ITS  1  OS  receiver  located  at  Gakona.  This  clock  uses  a  single¬ 
frequency  GPS  receiver  chip  set  that  can  also  provide  estimates  of  the  receiver  location.  We  use 
these  data  to  measure  the  position  error  in  these  estimates  by  comparing  the  instantaneous 
position  estimates  to  estimates  derived  from  long-term  averages  of  the  receiver  location.  The 
departures  of  the  instantaneous  positions  from  the  average  positions,  both  in  the  horizontal  and 
vertical  directions,  are  shown  in  the  plots.  These  plots  show  up  to  36  hours  of  data,  and  the 
horizontal  data  are  fitted  by  an  rms-based  error  ellipse. 

All  of  the  GPS  data  described  in  this  section  are  archived  at  either  the  HAARP  site  or  on 
servers  at  NWRA’s  Bellevue  office. 
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Figure  2.  Sample  updated  tomography  output  graphic.  The  top  panel  shows  the  tomography-derived 
electron  density  latitude-by-altitude  map  across  eastern  Alaska,  the  center  panel  the  vertical  TEC 
calculated  from  the  electron-density  profiles  output  by  the  tomography  processor,  and  the  lower  panel  is 
the  f0F2  estimates  from  the  profiles. 
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2.1.3.  Mobile  LEO  Satellite  Receiver  System 

The  NWRA  ITS30m  Receiver  System  was  conceived  as  a  mobile  system  wherein  the  user 
could  operate  from  an  ordinary  SUV-type  vehicle,  and  program  the  sample  rate.  The  concept 
was  to  modify  the  current  NWRA  ITS30  receiver  design  as  little  as  possible  to  allow  for  mobile 
operations.  The  system  is  also  designed  for  variable  data  collection  rates  to  allow  a  wider  range 
of  application,  particularly  for  observing  small-scale  ionospheric  irregularities. 

Variable  Sample-Rate  Design 

Current  series  of  ITS  receivers  have  a  fixed  sample  rate  of  50  Hz.  One  of  these  was 
modified  to  allow  for  collection  of  data  at  a  2,000-Hz  rate,  but  that  was  a  breadboard 
implementation  rather  than  a  design  modification.  It  is  simple  to  adjust  the  sample  rate  from  a 
computer,  as  the  A/D  card  used  to  digitize  the  input  analog  data  stream  is  under  computer 
control.  However,  as  in  any  data  sampling  scheme,  the  analog  signal  to  be  digitized  must  be 
filtered  to  prevent  aliasing.  In  the  case  of  a  50-Hz  sample  rate,  the  filter  is  a  25-Hz  4-pole  low- 
pass  filter.  The  design  effort,  then,  was  to  adjust  the  cutoff  frequency  of  the  anti-aliasing  filter  at 
the  same  time  as  the  digitizing  rate  was  changed. 

The  central  component  in  this  design  is  a  digitally  controlled  potentiometer  (pot).  The  filter 
parameters  are  set  with  capacitors  and  resistors,  and  so  a  digitally  controlled  resistor  can  be  used 
to  affect  the  necessary  tuning.  This  was  done  in  three  steps: 

1.  Identify  a  digital  pot  that  would  work  at  the  voltages  present  in  the  existing  circuit.  Many 
digital  pots  are  available,  but  the  great  majority  of  them  work  at  5  volts  whereas  the  existing 
circuit  works  at  +/-  15  volts.  One  such  pot  was  found  from  MAXIM,  an  electronic  component 
manufacturer. 

2.  Change  the  component  values  of  the  remaining  resistors  and  capacitors  in  the  filter  so  that  the 
filter  characteristics  (the  'shape'  of  the  filter  cutoff)  would  remain  the  same  at  any  value  of  the 
digital  pot.  This  was  accomplished  with  the  aid  of  a  circuit  simulation  program  from  Linear 
Technology,  a  leading  semiconductor  manufacturer. 

3.  Rework  the  circuit  board  artwork  to  accommodate  the  new  design. 

With  a  few  sample  boards  fabricated,  the  concept  was  tested,  and  a  number  of  design  defects 
were  discovered,  a  usual  occurrence.  The  test  board  was  modified  with  cuts  and  jumper  wires  to 
correct  the  defects  and  continue  testing  until  the  design  was  verified.  Based  on  the  testing,  new 
artwork  was  prepared. 

Note  that  this  design  was  specifically  adapted  to  the  use  of  a  certain  A/D  data  input  card 
from  National  Instruments,  the  6052E/6024E.  This  card  is  in  use  in  all  of  the  existing  ITS 
receiver  stations.  At  about  this  time,  we  began  to  consider  that  the  6052E/6024E  cards  had  been 
in  production  for  about  12  years,  which  is  about  as  long  as  hardware  of  this  kind  is  supported  by 
the  manufacturer.  Discussions  with  the  manufacturer,  National  Instruments,  indicated  that  the 
card  is  a  popular  one,  so  that  production  will  likely  continue  for  some  time  longer,  but  no  firm 
prediction  was  made  as  to  how  long  that  would  be.  In  view  of  the  above,  it  was  decided  to  take 
another  step  in  the  design  of  the  variable-sample-rate  system  to  incorporate  a  more  modem  A/D 
input  card,  the  USB  6009  (also  from  National  Instruments). 
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The  resulting  design,  now  in  production,  allows  the  Receiver  Detector  Module  (the  item  of 
receiver  hardware  that  is  the  focus  of  the  variable-sample-rate  concept)  to  be  built  to  interface 
with  either  the  6052E  or  the  USB  6009,  and  to  be  built  with  or  without  the  variable-sample-rate 
feature.  The  actual  receiver,  ITS30-15m,  is  fitted  with  a  Detector  configured  for  variable 
sampling  and  the  6024E. 

Mobile  Operation 

By  Mobile  Operation,  we  refer  to  the  operation  of  an  ITS  system  from  a  vehicle  such  as  an 
automobile  or  SUV.  To  achieve  mobile  operation,  three  issues  needed  to  be  addressed:  the 
computer,  power  and  the  antenna. 

1.  Mobile  Computer.  While  there  is  no  functional  reason  why  the  standard  ITS  configuration 
could  not  be  used  in  a  mobile  operation,  the  power  consumption  and  unwieldy  physical  aspects 
of  having  a  desktop  computer,  monitor,  keyboard  and  mouse  indicated  a  better  solution.  In  this 
we  were  fortunate  that  the  A/D  input  card  discussed  above  (6024E)  also  comes  in  a  PCMCIA 
version;  that  is,  a  version  that  fits  into  a  special  card  slot  that  is  standard  on  most  Laptop 
computers.  Therefore,  the  mobile  computer  is  simply  a  laptop  computer  with  ITS  software  and  a 
PCMCIA  6024E  installed. 

2.  Mobile  Power.  All  of  the  components  of  the  mobile  system,  receiver,  computer  and  GPS 
clock,  come  standard  operating  from  120  VAC.  In  another  stroke  of  good  fortune,  none  of  these 
components  are  overly  sensitive  to  the  'cleanliness'  of  that  power.  Therefore,  it  seemed  feasible 
to  create  a  power  system  using  the  host  vehicle's  12-volt  power  system,  and  a  12-volt  DC  to  120 
VAC  power  inverter.  Such  a  power  inverter  was  obtained,  and  the  system  was  tested  at  the 
Bellevue  laboratory  by  first  operating  from  120  VAC  from  an  extension  cord,  and  then  operating 
from  the  vehicle's  12-volt  DC  system  through  the  power  inverter.  No  difference  was  detected 
and  so  that  inverter  system  is  now  in  use. 

3.  Mobile  Antenna.  There  were  two  aspects  to  the  mobile  antenna  configuration:  the  adequacy 
of  having  only  the  vehicle  roof  as  a  groundplane  and  the  operational  issue  of  handling  the  ITS 
antenna.  Originally,  it  was  thought  that  a  large  ground-plane  structure  would  need  to  be 
fabricated  that  would  fit  atop  the  mobile  vehicle.  Considering  that  every  campaign  in  which  this 
system  was  likely  to  be  deployed  would  involve  a  different  (rented)  vehicle,  this  was  seen  as  a 
challenge.  It  was  then  observed  that  ITS  receivers  operating  in  Taiwan  were  consistently  using  a 
1-m  x  1-m  ground  plane.  This  gives  these  users  an  adequate  antenna  pattern,  and  gave  us  the 
idea  that  an  SUV  rooftop  might  be  sufficient.  Once  again,  testing  of  this  idea  indicated  that  at 
relatively  high  satellite  elevations,  the  SUV  roof  ground  plane  was  sufficient. 

Handling  of  the  ITS  antenna  system  is  awkward  because  the  antenna  mast  is  about  30  inches 
tall,  with  four  18-inch  elements  radiating  from  the  top  of  the  mast  structure,  four  7-inch  elements 
radiating  from  a  point  about  12  inches  from  the  base,  and  four  more  2-inch  elements  about  4 
inches  from  the  base.  The  solution  was  to  make  each  of  the  upper  8  elements  into  2  pieces 
threaded  together.  In  transport,  all  of  the  elements  project  only  2  inches  from  the  center  of  the 
mast.  When  set-up,  the  upper  and  middle  element  sets  have  appropriately  sized  element 
extensions  screwed  into  the  2-inch  stubs  to  form  a  full-sized  antenna.  The  mast  itself  is  affixed 
to  an  18”  x  24”  flat  aluminum  plate  which,  in  turn,  is  placed  on  the  vehicle  rooftop. 
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Testing  and  Operation 

Three  tests  of  the  mobile  system  were  conducted.  The  first  two  were  developmental  testing, 
and  were  not  intended  to  have  any  scientific  value.  The  final  test  involved  mobile  operations 
that  were  integrated  into  the  science  campaign. 

The  first  test  was  conducted  at  the  HA  ARP  site  at  a  location  approximately  50  feet  away 
from  the  fixed  ITS  station  there.  It  demonstrated  that  the  concept  was  sound:  the  mobile  data 
was  roughly  similar  to  the  fixed  station  data.  It  also  revealed  a  circuit  design  defect  that  required 
further  testing  at  the  Bellevue  lab  to  locate  and  correct.  It  was  at  this  test  that  the  idea  of  having 
the  antenna  elements  be  removable  was  conceived. 

The  next  test  was  also  conducted  at  the  HAARP  site  environment,  this  time  in  actual  off-site 
deployment  mode.  The  receiver  performed  correctly,  but  two  more  operational  issues  were 
revealed.  The  first  was  that  synchronizing  the  time  on  the  mobile  computer  to  UTC  was  very 
awkward  and,  at  best,  inaccurate.  The  solution  was  to  include  a  GPS  clock  into  the  suite  of 
equipment  used  for  the  mobile  operation.  The  clock  used  was  the  CNS  clock  that  is  part  of  the 
Bellevue  lab  station,  and  a  small  GPS  patch  antenna. 

The  second  operational  issue  was  that  finding  the  remote  site  where  the  mobile  station  was 
to  be  located  was  very  clumsy  using  simply  a  GPS  position  indicator.  To  solve  this  problem,  it 
was  necessary  to  plot  the  satellite  path  on  a  map  that  also  showed  the  roadways  in  the  vicinity, 
and  then  actually  visit  the  remote  site  location  early  enough  to  find  a  suitable  place  to  pull  off  of 
the  highway  and  set  up  the  station.  By  the  last  test,  the  TOPO  software  package,  from  National 
Geographic,  was  used  which  allowed  the  satellite  ground  path  to  be  plotted  on  a  map,  and  then 
the  resulting  map  to  be  printed  out. 

The  final  test  was  a  fully  operational  part  of  a  science  campaign.  All  the  planned  sites  were 
located,  the  station  set  up  and  operated,  and  data  collected.  In  several  cases,  the  satellite  that  was 
the  subject  of  observation  was  either  not  transmitting,  or  had  an  inadequate  transmitting  antenna 
pattern.  Therefore,  not  all  the  data  sought  was  collected  or  useful.  Samples  of  data  collected 
from  this  mobile  receiver  during  the  campaign  will  be  presented  later  in  this  report. 

During  this  final  test,  another  operational  issue  was  discovered.  A  second  ITS30  system, 
owned  and  operated  by  personnel  from  AFRL,  was  deployed  alongside  the  NWRA  ITS30m 
system  during  the  campaign.  The  AFRL  system  was  different  in  that  it  was  not  a  variable  data 
rate  system,  and  it  made  use  of  a  different  antenna  system  than  the  NWRA  system.  Otherwise,  it 
had  a  similar  receiver  and  a  laptop  computer.  The  subject  operational  issue  was  discovered  when 
the  AFRL  people  placed  the  antenna  LNA/preselector  module  inside  the  vehicle.  It  is  normally 
positioned  close  to  the  antenna  but,  since  the  antenna  was  so  close  to  the  receiver  (in  this  mobile 
configuration),  the  module  was  brought  inside  the  SUV.  In  this  arrangement,  the  cables  from  the 
antenna  to  the  LNA  and  the  cables  from  the  LNA  to  the  receiver  were  close  together  and  formed 
a  feedback  loop.  No  signal  was  detectable  under  these  conditions.  (Note  that  there  were  other 
problems  with  the  AFRL  receiver,  so  that  system  would  not  have  performed  well  in  any  case,  but 
the  feedback  problem  was  instructive  to  the  future  operation  of  the  mobile  receiver  concept.) 

To  provide  assurance  that  the  feedback  situation  will  not  reoccur,  and  to  minimize  the 
feedback  signal  path,  the  cable  set  from  the  antenna  to  the  LNA  was  shortened  so  that  the  LNA 
can  only  be  operated  in  close  proximity  to  the  antenna. 
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Summary 

Development  of  the  ITS30m  system  involved  electronic  design  effort  that  met  the  initial 
goal,  variable  data  rate,  and  was  extended  to  provide  for  the  use  of  more  modem  vendor-supplied 
equipment  in  the  future.  Development  also  required  overcoming  several  nut-and-bolt  operational 
issues  to  create  a  functional  system.  Figure  3  shows  the  ITS30m  system  as  completed.  Figure 
3  a  shows  the  system  prior  to  deployment  (the  open  box  next  to  the  antenna  holds  the 
preamplifier  electronics),  3b  shows  the  inside  of  an  SUV  at  the  HAARP  site  during  the  February 
2008  campaign,  and  3c  shows  an  outside  view  of  the  SUV  illustrating  the  antenna  and 
preamplifier  box,  also  from  the  February  2008  campaign. 


a.  b.  c. 


Figure  3.  The  ITS30m  system,  (a)  The  entire  ITS30m  system  prior  to  deployment,  (b)  The  system  deployed 

in  an  SUV.  (c)  The  antenna  and  preamplifier  box  as  deployed  for  pass  acquisition. 

2.2.  Tomographic  Imaging  of  the  Auroral  E  Region 

The  tomography  analysis  system  used  to  invert  observations  of  relative  TEC  from  the 
HAARP  receiver  chain  to  obtain  altitude -by-latitude  maps  of  ionospheric  electron  density  was 
developed  with  a  focus  on  imaging  the  major  features  of  the  F2  and  topside  ionosphere.  While 
there  was  anecdotal  evidence  that  the  process  was  capable  of  producing  images  that  showed 
realistic  structures  in  the  E  region,  specifically  auroral  E,  we  had  not  taken  a  systematic  look  at 
any  specific  cases  to  see  if  the  results  were  consistent  with  other  observations.  In  preparing  for 
the  February  2008  campaign  (discussed  elsewhere  in  this  report),  we  found  two  tomography 
images  from  passes  roughly  90  minutes  apart  that  showed  a  very  strong  E  region  in  the  initial 
image  with  a  much  weaker  E  region  in  the  second.  We  have  taken  a  closer  look  at  these  two 
images  and  compared  them  with  data  collected  from  other  HAARP  diagnostic  instruments 
during  the  times  of  the  two  passes. 

Figures  4  and  5  show  the  images  generated  from  the  data  collected  by  the  HAARP 
tomography  chain.  Figure  4  is  an  image  generated  from  an  Oscar  32  pass  over  the  time  interval 
1226-1246  UT  on  10  February  2008,  and  Figure  5  is  an  image  from  an  Oscar  25  pass  over  the 
time  interval  1358-1414  UT  on  10  February  2008.  An  intense  E-layer  stretching  from  Gakona 
(Gak)  to  Arctic  Village  (ArV)  is  very  evident  in  Figure  4,  with  a  very  weak  E-layer  over  the 
southern  end  of  this  latitude  range  in  Figure  5.  The  working  hypothesis  is  that  the  E-layer  in 
Figure  4  is  due  to  very  intense  auroral  precipitation  which  has  weakened  markedly  by  the  time  of 
the  second  pass  shown  in  Figure  5. 

The  images  in  Figure  6  are  a  series  of  Digisonde  soundings  take  from  the  Gakona  (left 
column)  and  College  (right  column)  Digisondes  during  the  times  covered  by  the  tomography 
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images  shown  in  Figures  4  and  5.  The  soundings  start  at  1200  UT  on  10  February  and  run  in  15 
minute  steps  through  1415  UT,  skipping  the  sounding  at  1315  UT  (the  Gakona  1315  UT 
sounding  is  nearly  identical  to  the  Gakona  1330  UT  sounding,  and  the  College  1315  UT 
sounding  is  nearly  identical  to  the  College  1300  UT  sounding).  The  soundings  prior  to  1230  UT 
show  very  substantial  E  layers,  followed  by  three  successive  soundings  where  absorption  is  so 
high  that  no  discemable  echo  trace  can  be  found  in  either  sounding.  This  is  during  the  period  of 
the  image  in  Figure  4,  which  shows  very  strong  E  layers  over  both  Digisondes.  The  E-layer 
traces  are  beginning  to  return,  and  by  the  time  of  the  second  pass  (Figure  5),  the  E  layer  over 
Gakona  has  weakened  and  there  is  an  F2-layer  echo  beginning  to  show  at  a  critical  frequency  of 
around  3  MHz,  which  is  consistent  with  the  densities  shown  in  Figure  5. 

This  picture  of  an  intense  auroral  precipitation  event  is  supported  by  the  data  shown  in 
Figures  7  and  8,  which  are  24-hour  plots  of  the  absorption  measured  by  the  HAARP  30-MHz 
classic  riometer  and  the  H,  D,  and  Z  traces  from  the  HAARP  magnetometer.  The  riometer 
shows  very  strong  absorption  at  the  time  of  the  first  pass  (Figure  4),  at  the  same  time  the 
magnetometer  is  showing  a  large  negative  bay  in  the  H  and  Z  traces.  The  riometer  absorption 
has  decreased  significantly  by  the  time  of  the  second  pass  (Figure  5),  as  has  the  intensity  of 
auroral  currents  as  observed  in  the  magnetometer  data.  Unfortunately,  there  were  no  DMSP 
passes  near  the  HAARP  facility  at  the  time  of  this  event,  but  data  from  other  longitude  sectors 
showed  strong  auroral-zone  precipitation  around  the  time  of  the  first  pass. 

We  plan  to  pursue  this  ability  to  image  the  auroral  E  region  with  the  tomography  system  in 
our  future  work  on  the  new  contract. 


9 


Altitude  (km) 


Geomagnetic  Latitude  (deg),  defined  at  350-km  altitude 

Figure  4.  Tomography  image,  Oscar  32  pass  from  1226  to  1246  UT  on  10  February  2008. 


1000 


800 


£  600 

<u 

a 

|  400 


200 


0 

56  59  62  65  68  71  74  77 

Geomagnetic  Latitude  (deg),  defined  at  350-km  altitude 

Figure  5.  Tomography  image,  Oscar  25  pass  from  1358  to  1414  UT  on  10  February  2008. 


2008Febl0;  1358  - 1415  UT;  North  to  South 


Cor  Gak  Del  ArV 

i  A  |  A  A  .  A 

j _ i _ l _ i _ i _ l _ i _ i _ 1 1 _ i _ i _ l _ i 


10 


\  Lowell 
'  DIGISONDE 


EoF2 

EoFl 

EoFlp 


foEs 


N/A 

N/A 

N/A 

N/A 

.35 

N/A 

5.03 


KUF(D) 

M(D) 


N/A 

N/A 


h'F 
h  F2 


N/A 

N/A 


N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 


Auto: 

Artist5 

500200 


6A76;_;008041i:00 

f  Lowell 
'  DIGISONDE 


»  D  :i5-0  E  Ion-Png  v. 

[  Pi  FFS  S  AX1I  PPS  IGA  PS 
1  713  100  20+  A1 


£oF2 

EoFl 

EoFlp 


N/A 

N/A 

N/A 

N/A 

.35 

N/A 

N/A 


MUF(D) 

M(D) 


N/A 

N/A 


h'F 

h'F2 


N/A 

N/A 


hmF2 

hmFl 


N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 


Auto: 

Artist5 

500200 


&*76;_;oo804ii;is 

\  Lowell 
'  DIGISONDE 


L'K'LV 

* 

' 

' 

:  jl 

F  . ^ 

foF2 

foFl 

EoFlp 


EoEs 


N/A 

N/A 

N/A 

N/A 

.35 

N/A 

N/A 


MUF(D) 

M(D) 


N/A 

N/A 


h'F 
h  F2 


N/A 

N/A 


N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 


Auto: 

Artist5 

500200 


ITHM 

L’L'IV 

f  Lowell 
’  DIGISONDE 

foF2 

EoFl 

EoFlp 


yF2 

yFl 


D  100  200  400  600 

HOF  .0  .0  .0  .0 

C0764_200804112000S.MMM  /  150£ 

f  Lowell 
’  DIGISONDE 


EoF2 

EoFl 

EoFlp 


N/A 


700 


HOF(D)  N/A 

H(D)  N/A  500 

D  3000.0 


•:::  ’  • 

L- 

.  ."Ill:-  :: 

'  1 

;  ::: 

•  . 

mm 

.  pfijSpfe 

i  rn*  i  ::*r  > 

MW  i| 

Issi!!  v:"';  1 

...  ....  . 

D  100  200  400  600 

HOF  .0  .0  .0  .0 

C0764_2008041121505 .MHM  /  150£ 

f  Lowell 
'  DIGISONDE 


foF2 

EoFl 

foFlp 


N/A 


yF2 

yFl 


700  + 

650  : 

600  : 

550  : 

500- 

. __1_  .1:  J  .  \  -l...  J_1_L  •  .........  ~L... 

ESSE 

Aif 

H 

450  : 

. |  r  frr  j........-....:.:...-  jy‘  -  |— — £ | . -j— ViT 

400^ 

350: 

. I |  . .Jl — . 

300: 

.  r  4  4 !•  •.:•••=[ 

250: 

_ j_ . . . . !•  ••'; _ il _ i  | . '  '•  | . . _ [  j 

200: 

|  f  44  rf'' 

150: 

fr :  -  !•/  7r  4 •••••  -j  -4  •  -J . 

100: 

. . . . 4r .  :  ;;  .  ■  \  ■  -i  fr 

60- 

- 1 - 1 - 1 - — : - 1 - 1 - 1 - 

76;_;oo804ii:30 


D  100  200  400  600 

HOF  .0  .0  .0  .0 

C0764_2008041123005.MMM  /  1502 


Figure  6.  Digisonde  soundings  from  Gakona  (left  column)  and  College  (right  column). 
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Figure  6.  (continued). 
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2.3.  Tomographic  Image  Uncertainty  Specification 

2.3.1.  Background 

In  Section  2.2.  of  SRI  for  this  contract,  we  described  our  approach  to  determining  the 
uncertainty  in  density  and  f0F2  derived  from  the  tomographic  inversion  process.  That  approach 
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began  with  the  theoretical  specification  of  uncertainty  provided  by  the  framework  we  have  used 
for  the  inversion  process  at  the  core  of  the  tomographic  system,  and  expanded  on  that  by 
comparison  with  a  few  observations  of  f0F2  made  concurrent  with  our  tomography  images  from 
the  Digisonde  located  at  Gakona.  During  the  second  part  of  the  contract,  we  expanded  this  to 
include  an  assessment  of  the  uncertainty  in  the  vertical  TEC  derived  from  our  tomographic 
images,  and  ran  a  more  detailed  comparison  with  a  much  larger  set  of  Digisonde  soundings. 

2.3.2.  Assessment  of  Uncertainty  in  Vertical  TEC 

A  different  approach  was  used  in  assessing  the  uncertainty  in  vertical  TEC  derived  from  the 
ITS  images.  In  previous  status  reports,  we  described  how  observations  of  f0F2  from  the  Gakona 
Digisonde  were  used  to  assess  the  accuracy  of  the  theoretical  estimates  of  the  uncertainty  in  the 
images,  and  we  discovered  that  the  theoretical  estimates  were  both  too  low  and  did  not  account 
for  an  offset  in  the  uncertainty  that  reflects  the  tendency  of  the  analysis  to  provide  estimates  of 
f0F2  that  are  systematically  low.  There  are  no  corresponding  independent  observations  of  TEC 
available  for  a  similar  comparison.  Two  potential  sources  were  (1)  TEC  calculated  from 
Digisonde  measurements  of  the  F2  peak,  and  (2)  TEC  from  the  Gakona  GPS  receiver  system. 
Both  were  rejected  due  to  uncertainties  in  both  quantities  that  were  at  least  on  order  of,  if  not 
larger  than,  the  expected  uncertainty  in  the  ITS-derived  TEC.  The  Digisonde-derived  TEC  was 
based  on  an  estimate  of  the  topside  content  (between  hmF2  and  the  nominal  satellite  altitude  of 
1,000  km),  thus  potentially  over-  or  under-estimating  the  TEC  in  ways  that  are  difficult  to 
quantify.  The  GPS  observations,  on  the  other  hand,  include  content  from  1,000-km  altitude  out 
to  the  GPS  satellites  at  (nominally)  20,000-km  altitude.  While  the  TEC  in  this  height  range  is 
small,  it  can  be  on  the  order  of  2-5  TEC  units  (1016  el/m2).  Additionally,  there  is  inherent 
uncertainty  in  the  TEC  measurements  due  to  differential  group  delay  in  the  satellite  and  receiver 
which,  while  it  is  calibrated  using  the  SCORE  process,  has  uncertainties  on  the  order  of  a  few 
TEC  units  and,  at  the  high  latitudes  sampled  from  Gakona,  can  be  larger. 

Rather  than  deal  with  these  additional  uncertainties,  it  was  decided  to  use  simulated 
ionospheric  images  as  a  basis  for  this  study.  These  images  were  generated  using  the  simple  three 
Chapman  layer  profde  model  known  as  RBTEC  which  was  developed  in  the  mid-1970s  by 
personnel  at  the  Air  Force  Global  Weather  Central  ( Damon  and  Hartranft,  1970;  Flattery  and 
Ramsay,  1975).  Three  different  types  of  images  were  generated:  a  smooth  ionosphere  with 
either  no  latitude  gradients  or  with  a  constant  south-to-north  gradient  in  f0F2;  and  two 
modifications  of  these  smooth  cases  with  a  Gaussian-shaped  enhancement  in  f0F2  located  within 
the  image  range:  one  with  the  enhancement  situated  near  Gakona  and  the  other  with  the 
enhancement  situated  near  the  northern  end  of  the  chain.  Figure  9  shows  an  example  of  one  of 
these  “blob”  images.  TEC  was  calculated  for  the  six  stations  in  the  HAARP  chain  by  integrating 
through  these  images  to  pre-calculated  satellite  positions,  and  relative  TEC  was  produced  by 
subtracting  off  the  TEC  value  of  the  first  point  in  the  simulated  pass.  These  relative  TEC  records 
were  processed  using  the  standard  ITS  processing  with  each  of  the  three  EOF  sets  (day,  dusk, 
and  night)  currently  used  in  generating  images  at  HAARP.  Figure  10  shows  the  resulting  ITS 
reconstruction  of  the  image  shown  in  Figure  9  using  the  “Day”  EOF  set.  (Note  that  the  barely 
discemable  uncertainty  bars  on  this  plot  were  generated  from  the  unmodified  ITS  variance 
estimates.  As  will  be  shown  below,  these  underestimate  the  true  expected  uncertainty.) 
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These  three  image  types  were  all  controlled  by  setting  four  parameters:  the  latitude  gradient, 
the  f0F2  value  at  the  southernmost  end  of  the  image,  hmF2  for  the  entire  image  (hmF2  was  not 
varied  in  latitude),  and  a  factor  that  controlled  the  topside  scale  height  [this  is  an  additional 
multiplicative  factor  added  to  the  topside  section  of  the  RBTEC  model  described  on  page  341  of 
Flattery  and  Ramsay  (1975)].  Table  1  summarizes  the  parameter  settings  used  in  generating  the 
576  images  generated  for  this  study.  Since  each  of  these  images  were  processed  using  each  of 
the  three  EOF  sets,  the  study  is  based  on  a  total  of  1,728  analyses  and  over  300,000  individual 
estimates  of  f0F2  and  TEC.  A  database  was  built  from  these  analyses  which  included,  for  each 
profile  in  each  image,  the  f0F2,  hmF2,  and  TEC  from  both  the  original  and  reconstructed  images, 
as  well  as  the  estimates  of  the  lower  and  upper  bounds  of  both  f0F2  and  TEC  as  derived  from  the 
ITS-calculated  variance  (denoted  in  this  report  by  -  and  +  superscripts  and  referred  to  as 
uncertainty  estimates). 

Table  1.  Values  of  parameters  used  in  generating  the  simulation  image  set.  The  Gradient  parameter  is  the 
percent  change  in  f0F2  across  the  image  from  the  initial  value  shown  to  the  value  in  the  final 
(northernmost)  profile.  The  topside  scale  factor  is  applied  to  factors  calculated  by  the  RBTEC  model. 


Parameter 

Values 

Gradient: 

None  ,  10%  ,  25%  ,  50% 

Initial  f0F2: 

2.5 , 4.5 , 6.5 , 9.5  MHz 

hmF2: 

225 , 275 , 325 , 375  km 

Topside  scale  factor: 

1.0, 1.3  ,  1.6 

The  two  plots  in  Figure  11  summarize  the  initial  findings  from  this  analysis.  These  are 
histograms  of  the  percent  difference  between  the  f0F2  and  TEC  calculated  from  the  original 
images  and  those  calculated  from  the  ITS  reconstructions.  These  histograms  include  the  results 
from  all  1,728  analyses.  The  solid  curves  are  the  calculated  histograms,  the  dotted  curves  are 
normal  distributions  fit  (in  a  least-squares  sense)  to  the  histograms,  and  the  vertical  long-dash 
lines  indicate  the  averages  of  the  difference  between  the  original-image  values  and  the  X'  and  X+ 
values  (where  X  is  f0F2  or  TEC).  These  provide  estimates  of  the  expected  average  square-root 
of  variance  based  on  the  ITS  calculation  of  the  variance  for  the  image.  The  horizontal  long-dash 
line  indicates  the  “i/e”  point  in  the  histograms,  which  provides  a  location  on  each  histogram  that 
should  correspond  to  the  square-root  of  variance  (assuming  that  the  histogram  is  well  represented 
by  a  normal  distribution).  If  this  horizontal  line  crosses  the  histograms  at  the  same  locations  as 
the  vertical  horizontal  lines,  this  would  indicate  that  the  observed  variance  (for  these  300,000+ 
cases)  was  well  represented  by  the  modeled  variance.  As  can  be  clearly  seen  in  Figure  11,  this  is 
not  the  case.  Not  only  is  the  variance  too  low,  but  it  does  not  capture  the  offsets  from  zero 
clearly  seen  in  both  histograms.  This  is  consistent  with  what  was  discovered  in  comparing  with 
the  Gakona  Digisonde  observations. 

During  the  early  stages  of  this  analysis,  the  ITS  processor  was  run  using  the  modified 
variance  calculation  based  on  the  Gakona  f0F2  analysis  reported  earlier.  We  encountered 
difficulties,  however,  in  obtaining  consistent  results  with  the  uncertainties  in  f0F2  and  NmF2 
because  the  modifications  to  the  variance  used  in  generating  the  uncertainty  estimates  were  made 
partially  in  terms  of  plasma  density  and  partially  in  terms  of  plasma  frequency.  For  the  purpose 
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e  9.  One  of  the  simulation  electron  density  images  generated  for  the  uncertainty  study.  The  upper 
3anel  shows  a  contour  plot  of  the  electron  density,  the  middle  panel  is  the  vertical  TEC  generated  by 
ntegrating  through  the  image  in  the  upper  panel,  and  the  lower  panel  is  f0F2  derived  from  the  profiles 
n  the  image. 


Figure  10.  Results  from  the  ITS  reconstruction  of  the  image  shown  in  Figure  9.  In  the  middle  and  lower 
panels,  the  TEC  or  f0F2  from  the  original  image  is  shown  as  a  dotted  curve  and  that  calculated  from  the 
ITS  reconstruction  is  shown  as  a  solid  curve.  The  error  bars  are  uncertainties  in  TEC  and  f0F2 
calculated  from  the  theoretical  variances. 
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Histogram  of  dfoF2  -  All  Cases  (Orig  -/+) 


Histogram  of  dTEC  -  All  Cases  (Orig  -/+) 


Figure  11.  Histograms  of  the  difference  between  the  original  and  ITS  reconstruction  foF2  (left  panel)  and 

TEC  (right  panel).  See  text  for  description  of  the  long-dash  lines. 

of  this  study,  we  disabled  the  adjustments  we  had  implemented  as  a  result  of  our  comparison 
with  the  Gakona  f0F2  observations,  which  is  why  the  uncertainty  estimates  reflected  in  the 
location  of  the  vertical  long-dash  lines  in  Figure  1 1  for  f0F2  are  smaller  than  will  account  for  the 
variance  shown  in  the  histograms,  and  are  symmetric  about  zero. 

The  adjustments  made  to  the  ITS-calculated  variances  in  order  to  bring  them  more  in  line 
with  the  variances  consistent  with  the  histograms  shown  in  Figure  1 1  were  made  following  the 
same  general  algorithmic  approach  as  that  used  in  adjusting  the  variances  in  the  Gakona  f0F2 
study,  with  some  small  differences.  As  in  the  Gakona  study,  the  variance  calculated  from  the 
ITS  theory  (in  program  stat_comp)  is  multiplied  by  a  factor  of  the  form 

M=  (1  +  aN^Fl)  [1] 

where  NmF2  is  the  image-averaged  value  of  the  density  at  the  F2  peak  and  a  is  a  free  model 

parameter.  This  factor  is  applied  to  all  variances  in  the  variance  image  generated  by  the  analysis. 
The  offsets  in  the  uncertainty  values  (obviously  required,  as  can  be  seen  in  the  offset  and 
asymmetric  histograms  in  Figure  11)  are  then  applied  when  the  estimates  of  profile  parameters 
are  generated  (program  stat2image).  The  offset  calculation  is  implemented  as  a  second  free 
parameter,  b,  that  is  multiplied  times  either  NmF2  (for  the  f0F2  uncertainties)  or  TEC  to  provide 
the  necessary  offset  in  the  negative  and  positive  uncertainty  values.  The  desire  is  to  find  a  value 
for  the  a  parameter  in  this  equation  which  provides  good  estimates  of  the  uncertainty  of  all 
parameters  shown  in  Figure  11.  It  is  likely  that  the  offset  parameter,  b,  will  be  different  for  the 
peak  and  TEC  parameters  as  evidenced  by  the  fact  that  while  the  peak  parameters  show  that  the 
ITS  analysis  typically  underestimates  the  values  of  the  peak  parameters,  it  typically 
overestimates  TEC.  We’ll  identify  the  offsets  as  bp  and  bi  for  the  peak  and  TEC  calculations, 
respectively. 

It  was  discovered  that  while  we  were  able  to  match  the  uncertainties  as  embodied  in  the 
histograms  shown  in  Figure  1 1  for  f0F2  and  NmF2  using  the  same  values  for  a  and  bP,  this  still 
left  the  uncertainty  calculated  for  TEC  too  low.  To  account  for  this,  we  introduced  another 


19 


parameter,  a/,  which  is  a  simple  multiplicative  factor  applied  to  the  uncertainty  calculated  from 
the  ITS  variance  image  after  it  has  been  modified  using  the  M  factor  described  above.  After  a 
second  iterative  process,  we  found  that  the  TEC  uncertainty  calculated  from  the  variance  image 
generated  using  the  M  factor  needed  to  be  increased  by  roughly  80%  to  achieve  the  variance 
shown  in  the  histogram.  The  initial  values  derived  for  the  four  free  parameters  were: 

a  =  0.8  ,  a/=1.8  ,  6p=-0.33  ,  ^=+0.03 

The  results  of  these  modifications  on  the  uncertainties  calculated  are  shown  in  Figure  12,  which 
is  a  repeat  of  the  histograms  from  Figure  11  with  the  new  average  uncertainties  (the  vertical 
long-dash  lines)  shown.  Note  that  these  now  intersect  the  histograms  at  the  1/e  points  as  desired. 

When  these  parameters  were  used  to  generate  new  error  bars,  we  discovered  that  the 
calculated  error  bars  did  not  appear  to  be  representative  of  some  of  the  simulated-image  results. 
In  particular,  there  appeared  to  be  problems  with  some  of  the  reconstructions  made  using  the  Day 
EOF  set.  The  histograms  shown  in  Figures  11  and  12  were  regenerated  independently  for  the 
reconstructions  from  each  of  the  three  EOF  sets,  with  the  results  shown  in  Figure  13.  We 
decided  that  rather  than  fitting  to  the  aggregate  histograms,  we  would  set  the  lower  and  upper 
ranges  of  the  uncertainty  to  fit  the  lowest  and  highest  1/e  points  of  the  three  histograms.  This 
altered  the  values  of  the  four  parameters  to  the  following: 

a  =  1.3  ,  a,  =1.5  ,  bp=-0A0  ,  ^=+0.03 

This  set  provided  better  overall  representation  of  the  uncertainty,  and  the  value  for  a  is  now  very 
close  to  that  determined  in  the  analysis  of  the  Gakona  f0F2  Digisonde  data  (1.2).  This  also 
reduced  the  amount  by  which  the  TEC  uncertainty  is  increased  over  the  increase  needed  to  fit 
f0F2  from  80%  to  50%. 


Histogram  of  dfoF2  -  All  Cases  (New  -/+)  Histogram  of  dTEC  -  All  Cases  (New  -/+) 


Figure  12.  Histograms  as  shown  in  Figure  11  with  uncertainty  locations  from  the  analysis  of  the  aggregate 
histograms. 

The  need  to  increase  the  variance  by  a  factor  of  50%  for  the  TEC  uncertainties  to  match  the 
histograms  generated  from  the  simulation  studies  arises  from  at  least  two  sources:  (1)  a 
breakdown  in  an  underlying  assumption  made  in  calculating  the  TEC  variance  from  the  density 
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variance  map,  and  (2)  a  breakdown  in  the  implicit  assumption  that  the  multiplicative  factor 
which  correctly  adjusts  the  variance  at  the  F2  peak  (the  M  factor)  is  valid  for  all  heights  in  the 
profile.  We’ll  address  these  separately. 


Histogram  of  dfoF2  -  All  Cases  (New  -/+)  Histogram  of  dTEC  -  All  Cases  (New  -/+) 


(ofoF2-ifoF2)/ifoF2  (%)  (oTEC-iTEC)/iTEC  (%) 


Figure  13.  Histograms  as  shown  in  Figures  11  and  12  with  uncertainty  locations  from  the  analysis  of  the  Day, 

Night,  and  Dusk  EOF  set  histograms.  The  aggregate  histogram  is  the  solid  curve,  and  the  individual 

EOF  sets  are  dotted  curves. 

The  expected  uncertainty  in  vertical  TEC  derived  from  the  ITS  images  is  calculated  from  the 
density  variance  images  produced  by  the  ITS  software  by  summing  up  the  product  of  the  size  of 
the  vertical  bins  squared  times  the  variance  in  each  bin  of  the  vertical  profile.  This  provides  an 
estimate  of  the  total  variance  to  be  expected  in  the  TEC  based  on  the  density  variances  in  each 
profile  bin,  assuming  that  there  is  no  correlation  between  the  variations  of  Ne  in  each  bin.  The 
square  root  of  this  variance  is  then  used  to  provide  for  error  bars,  or  an  uncertainty  range,  around 
the  TEC  calculated  from  the  ITS  Ne  image.  An  underlying  assumption  used  in  deriving  this 
calculation  is  that  there  is  no  correlation  in  the  density  variations  in  the  various  bins.  If  there  is 
correlation,  extra  terms  show  up  in  the  equation  used  to  calculate  the  expected  TEC  variance 
from  the  density  variances.  Assessing  the  magnitude  of  this  source,  in  terms  of  additional  TEC 
variance,  is  difficult  using  the  simulation  set  we  have  employed.  However,  it  is  to  be  expected 
that  there  is  some  correlation  between  the  density  variations  in  adjacent  bins,  and  that  this  is,  at 
least  qualitatively,  a  clear  potential  source  for  “extra”  variance,  beyond  what  has  been  calculated. 

The  validity  of  the  assumption  that  the  M  factor  derived  at  the  F2  peak  can  be  applied  at  all 
altitudes  throughout  the  profde  is  clearly  called  into  question  by  the  results  shown  in  Figure  14. 
This  figure  shows  the  height  variation  of  the  average  square  root  of  variance  (a  good  estimate  of 
the  uncertainty)  for  a  randomly  selected  simulation  image.  The  heavy  solid  curve  is  that 
calculated  directly  from  the  initial  and  reconstructed  images,  the  light  solid  curve  is  calculated 
from  the  ITS-generated  variance  image,  and  the  dotted  curve  is  the  ITS-generated  variance 
scaled  to  match  the  variance  calculated  from  the  reconstructed  image  at  the  F2  peak.  The  two 
horizontal  lines  indicate  the  range  of  hmF2  in  the  reconstructed  image.  Two  things  are  clear  from 
this  plot:  (1)  the  height  variation  of  variance  as  generated  from  the  ITS  processor  is  not 
adequately  modeling  what  is  observed,  and  (2)  a  scaling  factor  that  brings  these  two  curves  into 
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agreement  at  the  F2  peak  is  not  valid  for  the  entire  profile,  particularly  just  below  the  F2  peak 
and  in  the  topside  of  the  profile,  and  is  likely  to  lead  to  an  underestimated  TEC  variance. 


0.0E+00  3.3E+04  6.7E+04  1.0E+05 

sqrt(variance)  el/m**3 

Figure  14.  Comparison  of  the  uncertainty  calculated  for  a  single  simulation  image  by  (1)  direct  comparison 
with  the  original  image  (heavy  solid  line)  and  (2)  calculated  from  the  ITS  theory  (thin  solid  curve).  The 
dotted  curve  is  the  theoretical  curve  scaled  to  match  the  direct-comparison  curve  at  the  F2  peak.  The 
horizontal  lines  indicate  the  range  of  the  height  of  the  F2  layer  in  the  image. 

Given  limitations  in  resources  available  to  continue  pursuing  this  aspect  of  the  ITS 
processor,  we  have  chosen  to  not  attempt  to  identify  which  of  these  two  sources  provides  how 
much  of  the  observed  underestimate  in  TEC  variance,  nor  to  search  for  other  potential  sources. 
We  will  implement  the  scheme  outlined  here  in  the  ITS  processor  and  use  it  to  calculate 
estimates  of  the  uncertainty  in  the  f0F2  and  TEC  values  which  are  posted  on  the  HAARP  web 
site.  Figure  15  shows  the  final  comparison  of  the  simulated  and  tomography  f0F2  and  TEC  using 
the  final  algorithm  for  the  uncertainty  bars.  Note  that  the  simulation  image  used  in  this  report  to 
illustrate  the  uncertainty  calculation  is,  despite  being  selected  at  random,  very  nearly  a  worst- 
case  example.  This  is  reflected  in  the  fact  that  the  differences  shown  in  Figure  15  between  the 
f0F2  and  TEC  reconstruction  and  the  “true”  values  from  the  simulated  image  are  at,  and  in  some 
cases  beyond,  the  uncertainty  bars. 
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2.3.3.  Comparisons  with  Digisonde  Observations 

In  the  original  investigation  of  uncertainty  in  our  estimates  of  f0F2,  a  small  sample  of  hand- 
analyzed  Digisonde  observations  of  f0F2  was  used  to  estimate  the  true  uncertainty  in  our  images. 
Subsequent  to  the  word  described  in  the  preceding  section,  we  were  able  to  obtain  a  much  larger 
set  of  post-processed  Digisonde  f0F2  and  hmF2  observations  from  both  the  Gakona  Digisonde 
and  a  second  Digisonde  located  at  College,  Alaska.  These  data  were  provided  by  Dr.  Ivan 
Galkin  of  the  Center  for  Atmospheric  Research  at  the  University  of  Massachusetts  Lowell 
campus.  Table  2  shows  months  for  which  data  were  provided  by  Dr.  Galkin  for  the  years  2001 
through  2006.  This  provided  over  13,000  observations  of  f0F2  over  this  six-year  period. 


Table  2.  Months  for  which  UML-supplied  Digisonde  data  are  available. 
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The  original  plan  for  these  data  was  to  revisit  the  algorithms  described  in  the  previous 
section  that  provide  estimates  of  uncertainty  in  f0F2  and  TEC  derived  from  our  tomography 
images  and  then  adjust  these  algorithms  as  indicated  by  the  analysis  results.  We  were  only  able 
to  begin  this  process,  as  with  the  completion  of  the  HAARP  facility,  the  focus  of  our  efforts 
shifted  to  support  for  campaign  operations  at  the  site  and  analysis  of  data  collected  during  those 
campaigns.  We  will  report  on  progress  made  to  date  on  this  effort.  It  will  serve  as  the  starting 
point  for  future  work  to  improve  our  estimation  of  uncertainty  in  products  derived  from  the 
tomography  images. 

Figures  16  through  21  summarize  the  results  of  the  comparison  as  conducted  to  date.  Each 
of  these  figures  consists  of  five  plots  showing  the  variation  of  the  difference  between  the 
tomography-derived  and  Digisonde-derived  f0F2  and  hmF2  as  a  function  of  day-of-year  (DOY) 
(Figures  16  and  17),  local-solar  time  (SLT)  at  the  tomography  array  center  longitude  (Figures  18 
and  19),  and  sunspot  number  (SSN)  (Figures  20  and  21).  Note  that  we  used  a  sunspot  number 
derived  from  the  observed  10.7-cm  solar  radio  flux,  which  is  a  better  surrogate  for  the  effect  of 
solar  EUV  radiation  on  the  ionosphere  than  is  the  observed  sunspot  number.  For  each  of  these 
six  figures,  the  results  have  been  divided  into  five  cases:  all  observations  (center  plot),  cases 
where  the  Day  EOF  sets  were  used  in  the  tomographic  analysis  (upper-left),  the  Night  EOF  set 
(upper-right),  dawn  hours  (local  midnight  to  noon)  using  the  Dawndusk  EOF  set  (lower-left), 
and  dusk  hours  (noon  to  midnight)  using  the  Dawndusk  EOF  set  (lower-right).  The  tomography 
system  automatically  determines  which  EOF  set  to  use  for  a  particular  satellite  pass  based  on  the 
solar  zenith  angle  at  the  center  of  the  tomography  array  at  the  center  of  the  pass  time.  Note  that 
the  gap  in  data  points  evident  in  the  SLT  plots  is  an  artifact  caused  by  how  times  are  assigned  to 
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the  tomography  results  and  are  not  due  to  a  gap  in  observations.  A  similar  thinning-out  of  data 
towards  the  latter  part  of  the  year  (as  seen  in  the  DOY  plots)  is  due  to  the  uneven  sampling 
provided  by  the  UML  Digisonde  data  set,  which  has  less  data  provided  from  the  months  at  the 
end  of  the  years  covered. 


Compare  foF2  (Both  ,  day) 


Compare  foF2  (Both  ,  night) 


Compare  foF2  (Both  ,  dawn) 


Doy 


Compare  foF2  (Both  ,  dusk) 


Figure  16.  Variation  of  the  difference  between  tomography-derived  and  Digisonde-derived  f0F2  as  a  function 
of  day  of  the  year  (DOY).  The  upper  two  plots  are  for  tomography  cases  in  which  the  Day  and  Night 
EOF  sets  were  used,  the  lower  two  plots  are  dawn  (lower-left)  and  dusk  (lower-right)  cases  using  the 
Dawndusk  EOF  set,  and  the  center  plot  is  all  cases.  The  heavy  line  is  an  average  calculated  in  15-day 
wide  bins,  and  the  two  light  lines  are  ±1  sigma  from  the  mean.  This  analysis  includes  data  from  both 
Gakona  and  College  Digisonde. 
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Figure  18.  Variation  of  the  difference  between  tomography-derived  and  Digisonde-derived  f0F2  as  a  function 
of  local-solar  time  at  the  center  of  the  tomography  array  (SLT).  The  figure  layout  and  plot  details  are  as 
in  Figure  16. 


27 


ITS  hmF2  -  Obs  hmF2  (km)  ITS  hmF2  -  Obs  hmF2  (km) 
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Compare  hmF2  (Both  ,  night) 
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Figure  19.  Variation  of  the  difference  between  tomography-derived  and  Digisonde-derived  hmF2  as  a 
function  of  local-solar  time  at  the  center  of  the  tomography  array  (SLT).  The  figure  layout  and  plot 
details  are  as  in  Figure  16. 
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Figure  20.  Variation  of  the  difference  between  tomography-derived  and  Digisonde-derived  f0F2  as  a  function 
of  sunspot  number  (SSN).  The  figure  layout  and  plot  details  are  as  in  Figure  16. 


TS  hmF2  -  Obs  hmF2  (km)  ITS  hmF2  -  Obs  hmF2  (km) 


Two  things  are  evident  from  even  a  cursory  look  at  these  results.  First,  we  should  adjust  our 
uncertainty  scheme  to  vary  with  the  EOF  set  employed  in  generating  the  tomography  image  and, 
perhaps,  add  in  variations  with  time,  season,  and  SSN.  Second,  we  should  revisit  our  algorithm 
for  determining  when  we  should  use  a  particular  EOF  set,  either  selected  from  the  three  we  have 
been  using  or  from  a  different  set  tailored  in  a  different  fashion.  Neither  of  these  activities  were 
possible  within  the  remaining  scope  of  this  contract,  but  we  will  lay  out  a  skeleton  plan  in  this 
report  for  addressing  these  in  the  future. 

The  results  shown  in  these  figures  can,  however,  provide  us  with  guidance  for  assessing 
uncertainty  in  F2-region  parameters  derived  from  the  tomography  images  using  the  present  set  of 
EOFs.  The  salient  points  to  note  are  as  follows: 

f0F2: 

1 .  As  was  seen  in  the  analysis  of  the  simulation  data  set,  the  tomography-derived  values  for 
f0F2  tend  to  be  systematically  lower  than  the  observed  values,  in  this  case  observed  by  a 
Digisonde.  However,  there  are  clearly  diurnal  and  seasonal  variations  in  this  offset,  which 
may  be  coupled,  and  also  the  behavior  is  different  for  the  different  EOF  sets. 

2.  The  performance  of  the  Dawndusk  EOF  set  is  clearly  much  better  on  the  “dawn”  side 
than  on  the  “dusk”  side;  performance  which  begins  to  worsen  as  daytime  conditions  are 
encountered  at  sunrise.  This  might  indicate  that  we  should  transition  sooner  to  the  Day  EOF 
set  on  the  dawn  side  and  later  from  Day  to  Dawndusk  on  the  dusk  side.  Also,  we  should 
look  at  using  a  different  EOF  set  for  evening  twilight  (dusk)  than  for  morning  twilight 
(dawn). 

3.  The  night  case  shows  the  most  scatter,  which  is  not  surprising  given  the  challenge  of  both 
imaging  the  night  auroral  ionosphere  and  of  measuring  the  F2  region  using  an  ionosonde. 

4.  There  is  a  clear  trend  in  the  average  curve  for  the  variation  with  SSN  for  the  Day  EOF 
case.  In  early  development  of  the  tomography  system,  particularly  in  design  and  selection 
of  the  EOF  sets  to  use,  we  speculated  on  the  need  to  tailor  the  EOF  sets  used  as  a  function  of 
the  point  in  the  solar  cycle.  The  results  shown  in  Figure  20,  as  well  as  the  corresponding 
result  for  hmF2  shown  in  Figure  21,  support  this  conjecture. 

hmF2: 

1.  There  is  clearly  a  diurnal  variation  in  the  performance  of  the  tomography  system’s  ability 
to  accurately  specify  the  height  of  the  F2  layer.  Both  the  Day  and  Dawndusk  EOF  sets  show 
a  tendency  to  overestimate  hmF2  during  the  day,  with  a  transition  to  a  nearly  zero  mean 
difference  between  the  tomography  and  Digisonde  values  at  night. 

2.  There  is  no  clear  seasonal  variation,  with  a  broader  spread  (in  terms  of  RMS  variation)  in 
the  difference  than  seen  in  the  SLT  plot. 

3.  While  the  scatter  in  the  Night  EOF  results  appears  larger  than  for  the  other  two  EOF  sets, 
this  is  not  as  large  a  difference  as  in  the  f0F2  results. 

4.  Similar  to  the  f0F2  SSN  variation,  the  Day  EOF  set  results  show  a  systematic  trend  that 
could  be  an  indication  that  we  should  investigate  different  Day  EOF  sets  as  a  function  of 
solar  cycle  epoch. 
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Figure  22  shows  histograms  of  the  f0F2  results;  in  this  case,  the  percent  difference  calculated 
from  (f0F2Digi  -  f0F2rrs)/foF2rrs.  Note  that  this  is  the  opposite  sense  of  the  data  shown  in  the 
scatter  plots,  which  show  (f0F2iTs  -  f0F2Digi).  This  change  is  made  to  facilitate  comparison  with 
the  simulation  results  discussed  in  the  previous  section.  In  particular,  the  upper-left  plot  in  this 
figure  (for  all  EOFs)  should  be  compared  to  Figure  13  in  that  section.  Two  things  are  evident  in 
the  results  shown  in  Figure  22.  First,  while  the  tomography  f0F2  values  tend  to  underestimate 
the  Digisonde  observations,  the  offset  is  much  less  than  that  shown  in  Figure  13.  Also,  the 
distribution  curve  is  a  bit  broader  than  that  obtained  from  the  simulation  data,  which  should  not 
be  surprising  given  that  the  simulation  result  was  obtained  from  a  smooth,  relatively  well- 
behaved,  model  ionosphere.  Assuming  we  continue  to  use  a  single  algorithm  for  all  EOF  sets, 
these  results  indicate  that  the  uncertainty  bars  for  f0F2  should  be  broader  than  are  being  produced 
by  our  current  algorithm  and  the  offset  of  the  error  bars  from  the  tomography  values  should  be 
less,  perhaps  much  less,  than  at  present.  A  quick  rerun  of  the  analysis  that  produced  the  values 
used  in  the  uncertainty  algorithm  described  in  the  previous  section  gave  values  of  3.4  for 
parameter  a  and  -0.20  for  parameter  bp,  reflecting  both  the  smaller  offset  (controlled  by 
parameter  bp )  and  the  broader  distribution  in  the  uncertainty  (controlled  mostly  by  parameter  a). 

An  interim  next  step,  based  on  retaining  the  current  simple  model  for  the  uncertainty,  is  to 
use  the  information  gleaned  from  this  preliminary  analysis  to  adjust  the  parameters  of  the  present 
model.  Since  the  model  for  uncertainty  in  the  TEC  derived  from  the  tomography  images  is 
coupled  to  the  model  for  f0F2  uncertainty  (through  the  adjustment  of  the  variance  output  by  the 
inverse  processor),  we  cannot  simply  implement  new  parameters  without  reconciling  the  changes 
to  the  TEC  uncertainty.  Since  we  have  no  source  of  independent  TEC,  other  than  GPS  TEC 
which  provides  TEC  over  a  much  larger  height  range,  we  will  need  to  make  this  reconciliation 
based  primarily  on  what  was  learned  from  the  simulation  analysis. 

The  steps  to  be  taken  to  make  use  of  the  comparison  of  the  tomography  and  Digisonde  data 
sets  to  improve  the  tomography  process  and  the  specification  of  uncertainty  in  the  resulting 
images  are  as  follows: 

1.  Make  a  more  careful  analysis  of  the  histograms  shown  in  Figure  22  to  derive  new  values 
for  the  f0F2  uncertainty  model  and  revisit  the  simulation  results  in  light  of  this  modification 
to  the  f0F2  uncertainty  to  reset  the  TEC  uncertainty  model  parameters. 

2.  Collect  additional  Digisonde  observations,  if  possible,  to  better  cover  the  last  several 
months  in  each  year  listed  in  Table  2. 

3.  Expand  the  analysis  described  in  this  section  to  better  separate  local  time  and  seasonal 
variations. 

4.  Review  the  use  of  EOF  sets  in  the  tomography  analysis.  This  will  include  refining  the 
algorithm  used  to  determine  which  set  to  use  (currently  tied  to  the  solar  zenith  angle)  and 
looking  at  gains  to  be  made  by  having  EOF  sets  tailored  by  SSN. 

5.  If  the  EOF  sets  are  changed,  or  their  use  is  changed,  reprocess  all  of  the  tomography  data 
sets  to  generate  new  images. 

6.  Redo  the  comparison  of  the  tomography-derived  f0F2  and  hmF2  with  the  Digisonde  data 
and  develop  a  new  model  for  uncertainty  based  on  this  new  analysis. 
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Percent  of  Obs  Percent  of  Obs 


Compare  Tomography  and  Digisonde  foF2 


Station:  both  EOF  Set:  All  Station:  both  EOF  Set:  day 


Station:  both  EOF  Set:  night  Station:  both  EOF  Set:  dusk 


Figure  22.  Histograms  of  the  percent  difference  between  Digisonde  f0F2  and  tomography  f0F2  for  all  cases 
(upper  left),  Day  EOF  set  (upper  right),  Night  EOF  set  (lower  left),  and  Duskdawn  EOF  set  (lower 
right).  Solid  curves  are  the  data,  and  the  dashed  curves  are  Gaussian  curves  fit  to  the  data. 


7.  Either  adjust  the  TEC  uncertainty  model  to  be  consistent  with  both  the  new  F2-region 
uncertainty  model  or  revisit  the  TEC  uncertainty  through  redoing  the  simulation  analysis 
with  the  new  EOF  sets. 
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The  goal  of  this  process  would  be  both  an  improved  performance  of  the  tomography 
processor  and  a  better  specification  of  the  uncertainty  in  parameters  of  interest  derived  from  the 
tomography  images. 

2.4.  Digisonde  Inclusion  in  Tomography 

As  we  have  previously  demonstrated,  the  tomographic  inversion  processor  we  have 
developed  is  capable  of  using  data  from  both  Digisonde  soundings  and  from  in-situ  density 
measurements  from  the  DMSP  SSIES  instrument  (see  4. 1.2. 3.2  in  Andreasen  et  al.,  2004). 
While  we  have  used  this  capability  in  post-analysis  inversions,  we  have  never  implemented  it  for 
use  in  our  near  real-time  inversions  which  are  posted  on  the  HAARP  web  site.  The  DMSP 
SSIES  data  are  not  available  until  well  after-the-fact  (we  obtain  these  data  from  Dr.  Fred  Rich, 
AFRL/VSBXP),  and  while  the  Gakona  Digisonde  data  are  available  in  near  real-time,  the 
difficulty  in  automated  analyses  of  high-latitude  ionospheric  soundings  and  lack  of  adequate 
quality-control  information  has  relegated  use  of  the  data  to  post-analyses  when  a  human  can 
review  the  Digisonde  analyses  for  accuracy.  We  have  found  that  using  either  type  of  additional 
data  as  input  can  markedly  improve  the  results  through  providing  “tie-down”  points  in  the 
vertical  profile,  not  just  in  the  vicinity  of  the  observation,  but  across  the  entire  image. 

We  have  begun  work  on  finding  a  way  to  be  able  to  use  at  least  some  of  the  Digisonde 
analysis  in  our  real-time  images.  As  shown  in  Figure  23,  there  is  gain  to  be  had  by  including  the 
Digisonde  f0F2  and  hmF2  estimates.  The  plots  on  the  left  side  of  Figure  23  show  the  tomography 
inversion  results  using  only  slant-range  relative  TEC  from  the  tomography  chain.  Note  the 
difference  between  the  f0F2  derived  from  the  tomography  inversion  (lower  panel  on  the  left  side) 
and  the  Gakona  and  College  Digisondes  as  indicated  by  asterisk  (*)  symbols.  This  tendency  of 
the  tomography  inversion  process  to  underestimate  f0F2  (and  overestimate  hmF2)  during  the 
daytime  hours  appears  to  be  fairly  systematic  (a  finding  of  our  ongoing  study  of  uncertainty  in 
the  tomographic  inversion  process).  The  plots  on  the  right  side  of  Figure  23  show  the  results 
when  the  f0F2  and  hmF2  from  a  single  Gakona  Digisonde  sounding  are  added  to  the  TEC  data. 
Note  that  the  input  of  this  single-point  observation  greatly  improved  the  f0F2  fit  at  the  College 
Digisonde  location,  which  was  not  input  to  the  analysis. 

Due  to  the  improvement  provided  by  using  the  Digisonde  data,  we  began  work  to  determine 
how  we  might  be  able  to  use  at  least  some  of  the  Digisonde  observations  in  our  real-time 
analyses.  Due  to  limitations  in  resources,  this  never  progressed  beyond  preliminary  discussions 
with  various  experts  in  the  field  of  quality  control  of  Digisonde  observations  (Dr.  Terry  Bullett, 
Mr.  Feo  McNamara,  Dr.  A.  Fee  Snyder,  Prof.  Bodo  Reinisch,  and  Dr.  Ivan  Galkin).  This  would 
be  a  fruitful  area  to  investigate  further. 

2.5.  Observations  of  the  Artificially  Perturbed  Ionosphere 

NWRA  was  involved  in  three  campaigns  at  the  HAARP  facility  near  Gakona,  AK,  during 
the  period  covered  by  this  report:  the  March  2006  HAARP  Optics  Campaign,  the  October  2008 
HAARP  Irregularities  and  Optics  Campaign,  and  the  February  2008  HAARP  Irregularities  and 
Optics  Campaign.  In  addition  to  running  experiments,  NWRA  was  given  the  responsibility  to 
coordinate  and  oversee  the  latter  two  of  these  campaigns.  The  October  2007  campaign  was 
scheduled  to  run  for  10  days  but,  due  to  difficulties  at  the  site,  this  campaign  was  terminated 
after  the  first  day.  No  results  were  obtained  from  that  campaign,  although  we  did  discover  an 
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Geomagnetic  Latitude  (deg),  defined  at  350-km  altitude 


Geomagnetic  Latitude  (deg),  defined  at  350-km  altitude 


Figure  23.  Tomography  results  for  12  November  2006.  The  results  in  the  three  panels  on  the  left  side  were 
generated  using  only  slant-range  relative  TEC,  while  the  results  shown  on  the  right  side  were  generated 
using  both  TEC  and  a  single  Gakona  Digisonde  observation. 

inconsistency  in  how  passes  were  identified  for  collection  between  the  fixed  and  mobile 
receivers  which  saved  us  from  losing  passes  in  the  February  campaign.  In  this  section,  we  will 
present  results  obtained  during  the  March  2006  and  February  2008  campaigns. 

2.5.1.  Results  from  March  2006  Campaign 

We  obtained  IRI  support  for  two  short  time  periods  during  the  March  2006  HAARP  Optics 
Campaign  in  conjunction  with  two  NIMS  passes  very  near  to  the  field-aligned  location  in  the 
heater  field-of-view.  A  major  difference  from  similar  experiments  attempted  during  the  2003 
optics  campaign  was  that  these  passes  occurred  during  morning  daylight  hours.  The  first  pass 
was  unsuccessful  due  to  sudden  onset  of  geophysical  activity,  generating  irregularities  from 
naturally  occurring  processes  and  producing  high  levels  of  absorption  in  the  lower  ionosphere, 
effectively  shielding  the  F2  region  from  the  heater  effects. 

We  were  more  successful  with  the  second  pass  on  27  March  2006,  with  the  time  of  nearest 
approach  to  the  field-aligned  point  at  00:33:51  UT.  Geophysical  conditions  were  very  quiet,  and 
we  were  able  to  observe  irregularity-generated  scintillation  on  the  VHF  and  UHF  signals  from 
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the  Oscar  25  satellite.  Data  from  the  SuperDARN  radar  located  at  Kodiak,  AK,  confirmed  that 
there  were  no  natural  irregularities  in  the  scale-sizes  the  radar  is  sensitive  to  both  prior  to,  and 
after,  the  time  period  during  which  the  heater  was  running.  Figure  24  shows  a  section  of  the 
VHF  and  UHF  detrended  signal  intensity  and  the  detrended  differential  phase  centered  on  the 
time  during  which  the  ray  path  passed  through  the  heated  region.  While  it  is  typical  to  see 
enhancements  in  scintillation  levels  when  the  raypath  approaches  the  field-aligned  point,  this 
enhancement  requires  the  presence  of  irregularities  to  operate.  Data  from  this  same  pass  at  Delta 
Junction  and  at  Cordova  both  showed  no  enhancement  in  scintillation  levels  at  their  respective 
field-aligned  points,  and  similar  passes  with  the  same  geometry  and  time-of-day  taken  at  Gakona 
have  also  not  shown  enhancements  of  this  level.  We  are  fairly  certain  that  the  irregularities  that 
caused  the  scintillation  shown  in  Figure  24  were  generated  by  the  heater.  Further  details  and 
analyses  of  this  pass  were  presented  in  a  poster  at  the  12th  Annual  Ionospheric  Interactions 
Workshop  in  Santa  Fe  (Secan  and  Fremouw,  2006). 


Gakona  I  27Mar06  0025UT  I  o25 


Figure  24.  Detrended  VHF  (top  panel),  UHF  (middle  panel),  and  differential  phase  (lower  panel)  from  a  200- 
second  segment  of  data  collected  at  Gakona,  AK,  from  the  Oscar  25  pass  at  0033  UT  on  27  March  2006. 
The  vertical  dashed  line  indicates  the  time  of  closest  approach  to  the  field-aligned  point  (as  observed 
from  Gakona),  and  the  horizontal  dashed  lines  indicate  the  time  during  which  the  ray  path  was  within  8° 
of  beam-center  (the  region  within  which  heater-generated  optical  emissions  are  typically  observed). 
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2.5.2.  Results  from  February  2008  Campaign 

As  noted  earlier,  NWRA  was  the  campaign  organizer  and  coordinator  for  this  campaign. 
Details  of  this  campaign,  including  a  discussion  of  lessons  learned  during  the  campaign 
operation,  have  been  provided  in  a  separate  report  ( Kennedy  et  al.,  2008).  In  this  report,  we  will 
focus  solely  on  the  results  obtained  from  NWRA’s  experiments  and  participation  in  experiments 
run  by  others. 

NWRA  was  assigned  seven  experiment  time  slots  during  this  campaign,  and  we  collaborated 
with  Dr.  Keith  Groves  of  AFRL  on  an  eighth.  Six  of  the  eight  were  LEO  satellite  passes  and  two 
were  joint  LEO  and  GPS  passes.  We  were  also  provided  time  originally  assigned  to  another 
experimenter  for  a  GPS  pass.  Table  3  provides  a  summary  of  the  experiment  times  NWRA  was 
provided.  In  this  table,  the  Case  column  specifies  the  case  number  we  assigned  the  experiment 
(the  Milikh  case  is  time  assigned  originally  to  Gennady  Milikh;  the  Groves  case  is  the 
experiment  we  ran  collaboratively  with  Dr.  Groves);  the  Satellites  column  indicates  which 
satellites  were  involved  (“G”  for  GPS  with  the  GPS  PRN,  “FM”  for  FORMOSAT,  “O”  for 
NIMS,  and  “F”  for  DMSP);  the  Program  column  indicates  the  configuration  of  the  IRI  beam;  the 
Freq  column  is  the  IRI  frequency;  and  the  UT  and  AST  columns  provide  the  date  and  start/end 
times  of  the  IRI  use.  The  three  programs  shown  were  STD  for  our  “standard”  program  of  CW 
O-mode  with  the  beam  directed  to  the  geomagnetic  zenith,  NS  for  a  beam  elongated  in  the  north- 
south  direction  aligned  with  the  array  geometry  (also  CW  O-mode),  and  30sOnOff  is  the 
standard  mode  with  a  30-seconds  on,  30-seconds  off  operation. 


Table  3.  NWRA  February  2008  Campaign  Experiment  Times. 


Case 

Satellites 

Program 

Freq 

MHz 

UT 

AST 

Date 

Start 

End 

Date 

Start 

End 

Milikh 

G21 

STD 

3.4 

21-Feb 

19:25:00 

19:45:00 

Thu, 21-Feb 

10:25:00 

10:45:00 

1 

FM6/G23 

STD 

2.84 

22-Feb 

03:47:00 

04:13:00 

Thu,  21-Feb 

18:47:00 

19:13:00 

2 

FM6 

STD 

4.5 

22-Feb 

20:22:00 

20:46:30 

Fri,  22-Feb 

11:22:00 

11:46:30 

3 

032 

STD 

5.2 

23-Feb 

22:40:00 

23:04:30 

13:40:00 

14:04:30 

4 

025 

NS 

4.2 

24-Feb 

01:31:00 

01:55:30 

Sat,  23-Feb 

16:31:00 

16:55:30 

5 

F15 

NS 

3.4 

24-Feb 

17:50:00 

18:15:00 

Sun,  24-Feb 

08:50:00 

09:15:00 

Groves 

023 

NS 

5.2 

24-Feb 

22:17:00 

22:47:00 

Sun,  24-Feb 

13:17:00 

13:47:00 

6 

FM5 

NS 

4.8 

25-Feb 

00:02:00 

00:25:00 

Sun,  24-Feb 

15:02:00 

15:25:00 

7 

F15 

STD 

2.7 

25-Feb 

03:44:00 

04:10:00 

Sun,  24-Feb 

18:44:00 

19:10:00 

7 

G23 

30sOnOff 

2.7 

25-Feb 

04:10:30 

04:18:00 

Sun,  24-Feb 

19:10:30 

19:18:00 

Data  were  collected  from  both  the  Gakona  fixed-asset  ITS  1  OS  receiver  and  the  NWRA 
ITS30m  mobile  receiver  for  most  of  these  time  periods.  Data  were  also  collected  from  NWRA’s 
HAARP  GPS  GIOS  receiver  for  Dr.  Gennady  Milikh,  and  have  provided  Dr.  Milikh  with  the  raw 
data  from  his  experiment  times. 

2. 5. 2.1.  Experiment  with  DMSP  F 1 5  Overflights 

Two  experiments  were  run  during  the  campaign  making  use  of  the  two-frequency  beacon 
and  in-situ  instrumentation  on  the  Defense  Meteorological  Satellite  Program  (DMSP)  FI 5 
satellite.  The  clearest  example  of  a  possible  feature  generated  by  the  HAARP  IRI  can  be  seen  in 
results  obtained  from  a  DMSP  F15  pass  near  the  HAARP  facility  around  0400  UT  on  25 
February  2008.  The  beacon  data  collected  from  this  pass  proved  to  be  unusable  for  extracting 
scintillation  information,  but  data  from  the  in-situ  instruments  shows  what  may  well  be  a  density 
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structure  at  the  satellite  altitude  (nominally  840  km)  generated  either  directly  or  as  a  by-product 
of  HA  ARP  IRI  heating  in  the  F-region. 

Figure  25  shows  the  geometry  of  this  pass  and  of  the  IRI  heated  region  as  configured  during 
the  time  of  the  pass.  The  left  panel  shows  the  track  of  the  DMSP  pass  as  mapped  down  field 
lines  (using  IGRF  2005)  to  an  altitude  of  220  km,  the  nominal  F2-layer  altitude  at  the  time,  and 
the  3-dB  and  6-dB  contours  of  the  IRI  beam  at  that  same  altitude.  The  heavy  line  indicates  the 
location  of  the  structure  we  will  be  studying.  The  right  panel  shows  the  satellite  track  at  the 
nominal  satellite  altitude  of  840  km  and  the  IRI  beam  contours  at  that  altitude.  Geomagnetic 
conditions  were  very  quiet  during  this  pass,  and  the  SuperDARN  radar  from  Kodiak  showed  that 
the  ionosphere  was  free  of  meter-scale  irregularities  in  the  vicinity  of  HAARP  aside  from  those 
generated  by  the  IRI. 
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Figure  25.  Track  of  DMSP  F15  field  line,  mapped  to  220-km  altitude,  for  pass  at  0400  UT  on  25  February 
2008  (dashed  and  solid  line).  Solid  section  of  this  track  indicates  the  section  of  patch  in  which  a  possible 
IRI-produced  feature  was  found.  The  two  ovoids  are  the  220-km  3-dB  and  6-dB  contours  of  the  IRI 
heater  as  configured  for  this  experiment. 

Data  from  the  F15  in-situ  instruments  are  shown  in  Figures  26a  (ion  density)  and  26b  (drift 
velocity  and  ion  temperature).  The  structure  of  interest  is  in  the  center  of  these  plots,  clearly 
seen  in  the  bottom  panel  of  Figure  26a.  This  is  an  isolated  structure  in  the  pass  data,  and  while  it 
does  not  map  directly  into  the  F-layer  region  in  which  the  IRI  is  interacting,  it  is  very  near  to  it. 
Our  tentative  hypothesis  is  that  the  larger-scale  structure,  as  seen  in  the  bottom  panel  of  Figure 
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26a,  is  plasma  that  has  been  moved  along  the  field  line  from  lower  altitudes,  with  the  finer-scale 
structure  evident  in  the  top  panel  of  Figure  26a  more  probably  generated  locally  due  to  plasma 
instabilities,  possibly  a  gradient-drift  instability.  The  fact  that  one  side  of  the  larger  structure  has 
this  finer-scale  structure  than  the  other  is  symptomatic  of  this  type  of  instability  in  isolated 
density  structures.  The  coincident  onset  of  upwelling  evident  in  the  vertical  drift  velocity  (center 
panel  of  Figure  26b)  is  consistent  with  this  picture,  as  is  the  decrease  in  ion  temperature  evident 
in  the  lower  panel  of  Figure  26b. 
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Figure  26.  Data  from  the  DMSP  F15  SSIES  instruments  from  the  pass  shown  in  Figure  25.  (a)  Data  from  the 
SSIES  Scintillation  Meter  (SM)  instrument.  The  bottom  panel  shows  the  total  ion  density,  the  top  panel 
the  output  of  a  digital  filter  applied  to  the  ion  density  data  with  a  band  pass  >  2.0  seconds  (a  detrended 
version  of  the  lower  panel),  and  the  middle  panel  shows  the  output  of  a  digital  filter  with  a  band  pass  of 
2.0  to  10.0  seconds,  (b)  Data  from  the  SSIES  Drift  Meter  (DM)  (top  two  panels)  and  Retarding 
Potential  Analyzer  (RPA)  (bottom  panel)  instruments.  The  top  two  panels  show  the  horizontal  and 
vertical  cross-track  ion  drifts,  respectively,  and  the  bottom  panel  shows  a  provisional  estimate  of  the  ion 
temperature  from  two  sources  (AFRL  and  NWRA  analyses).  In  all  panels,  the  heavy  vertical  line 
indicates  closest  approach  to  the  center  of  the  HAARP  IRI  heated  region,  the  heavy  horizontal  lines 
indicate  the  period  that  the  DMSP  satellite  was  within  the  3-dB  contours  of  the  heater  beam  (mapped  to 
840-km  altitude),  and  the  light  horizontal  lines  indicate  the  time  that  the  satellite  was  within  the  6-dB 
contours 

A  critical  question  that  must  be  answered  in  order  to  confirm  that  this  structure  is  indeed  due 
to  IRI  interactions  is  how  the  structure  was  formed.  Options  include  (1)  direct  interaction 
between  the  IRI  beam  and  the  ionosphere  either  at,  or  very  near,  the  satellite  altitude,  (2)  bulk 
transport  of  plasma  from  below  the  satellite  altitude  along  the  field  line  by  either  changing  the 
thermal  equilibrium  along  the  field  line  or  by  changing  the  electric  field  structure  along  the  field 
line,  or  (3)  some  form  of  “spoofing”  of  the  DMSP  instruments  by  direct  interaction  with  the  IRI 
RF  signal. 

While  all  of  these  have  strong  and  weak  points,  we  believe  that  option  3  is  unlikely  but 
should  be  verified  with  engineers  familiar  with  the  DMSP  instrumentation.  The  fact  that  the 
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horizontal  cross-track  velocity  shows  no  feature  where  the  vertical  velocity  shows  a  sharp 
change  in  the  flow  characteristics,  and  the  uneven  distribution  of  the  effect  on  all  parameters  in 
the  region  where  the  satellite  is  within  the  IRI  beam  at  840-km  altitude  (as  indicated  on  all  plots 
in  Figure  26),  are,  we  believe,  strong  circumstantial  evidence  that  this  is  unlikely  to  be  spoofing. 
Similarly,  we  also  believe  option  1  to  be  unlikely,  unless  a  plausible  mechanism  for  the 
interaction  can  be  proposed  given  the  IRI  operating  mode  and  frequency.  One  observation 
against  this  option  is  that  there  was  no  change  in  the  electron  temperature  as  measured  by  the 
Langmuir  probe  on  DMSP  (data  not  shown). 

We  believe  that  one  of  the  scenarios  in  option  2  is  the  most  likely,  although,  whether  the 
plasma  transport  is  due  to  thermal  or  electrodynamical  processes  is  not  clear  at  present.  If  the 
plasma  in  the  feature  has  been  transported  from  lower  in  the  ionosphere  without  any  loss,  if  we 
assume  a  scale  height  of  about  150  km  and  a  Chapman  profile  the  density  difference  between  the 
center  of  this  feature  and  the  density  away  from  the  feature  would  correspond  with  densities 
roughly  50  km  below  the  satellite  altitude  when  outside  the  feature.  With  an  upward  velocity  of 
200  m/s,  that  infers  a  transit  time  on  the  order  of  four  minutes,  which  is  well  within  the  time 
frame  of  the  experiment.  We  will  be  working  with  Dr.  Brenton  Watkins  and  his  student,  Chris 
Fallen,  at  the  University  of  Alaska  Geophysical  Institute  to  use  their  one-dimensional  model  of 
the  expected  behavior  of  the  dynamics  along  the  field  line  to  sort  out  the  different  processes  that 
could  lead  to  the  generation  of  this  structure. 

2. 5. 2. 2.  Dual-Receiver  LEO  Satellite  Observations 

NWRA  ran  a  series  of  experiments  with  a  variety  of  LEO  satellite  passes  using  both  the 
fixed-position  ITS  1  OS  receiver  located  at  the  HAARP  site  (designated  the  Gakona  receiver)  and 
the  new  mobile  ITS30m  receiver  (designated  the  Mobile  receiver)  which  was  deployed  in  a 
rented  SUV  at  a  number  of  locations  around  the  site.  AFRL  (Dr.  Keith  Groves)  also  had  a  new 
ITS30m  receiver  deployed  at  this  campaign,  which  was  deployed  in  tandem  with  the  NWRA 
Mobile  receiver  in  order  to  check  out  this  new  receiver  as  a  sort  of  shake-down  cruise.  The 
ITS30m  deployments  were  under  the  guidance  of  NWRA  consultant  Mr.  Frank  Smith,  who  also 
assisted  the  AFRL  team  that  deployed  Dr.  Groves’  ITS30m  receiver.  Mr.  Smith  also  spent  a 
considerable  amount  of  time  helping  the  AFRL  team  and  Dr.  Groves  sort  out  problems  with  their 
receiver  and  how  it  was  deployed.  A  problem  was  discovered  in  the  A/D  converter  card  in  the 
laptop  used  by  AFRL  for  data  collection  from  their  receiver,  and  other  issues  relating  to  antenna 
cabling  were  also  resolved. 

We  were  able  to  collect  useful  data  from  several  of  these  experiments,  one  of  which  will  be 
presented  later  in  this  report.  We  did  encounter  problems  with  several  of  the  satellites,  however, 
that  made  it  difficult  or  impossible  to  collect  useful  data  from  them: 

1.  We  discovered  that  not  all  FORMOSAT-3  (COSMIC)  satellites  have  their  beacons  turned  on 
when  these  satellites  pass  over  Alaska.  The  organization  that  controls  these  satellites  was 
requested  to  have  several  satellites’  beacons  enabled,  but  fewer  were  on  than  expected.  For 
future  campaigns,  the  status  of  these  beacons  needs  to  be  established  further  in  advance  of  the 
campaign  with  a  clearer  understanding  of  which  beacons  will  be  on. 

2.  The  signals  from  the  FORMOSAT-3  satellites  that  we  were  able  to  collect  showed  severe 
structure  due,  mostly  likely,  to  the  transmitter  antenna  pattern.  The  reduction  in  UHF  signal 
strength  in  a  deep  null  near  overhead  passes  by  these  satellites  caused  all  receivers  (the  fixed 
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ITS  1  OS  and  both  mobile  receivers)  to  lose  lock  on  the  satellite  as  it  reached  overhead.  Even  on 
passes  where  the  signals  were  tracked  for  the  entire  pass,  removing  the  artificial  structure  from 
the  antenna  pattern  from  the  received  signals  in  order  to  study  the  IRI  heater-produced  structures 
will  be  challenging.  For  this  reason,  we  plan  to  place  a  lower  priority  on  FORMOSAT-3  passes 
in  future  campaigns. 

3.  A  similar  problem  with  structure  was  found  in  the  DMSP  F15  beacon  data.  We  suspect  that 
this  is  a  combination  of  both  structure  in  the  transmit  antenna  pattern  and  interference  with  the 
signal  by  structures  on  the  bottom  of  the  DMSP  spacecraft  that  can  block  or  reflect  the  beacon 
signal.  We  also  found  small-scale  oscillations  (order  1-2  seconds)  in  both  the  phase  and  intensity 
from  this  beacon,  which  will  make  the  data  very  difficult  to  use  for  scintillation  observations. 
These  issues  have  been  identified  to  Dr.  Paul  Bernhardt,  NRL,  who  is  the  PI  for  this  beacon. 

One  of  the  issues  we  wanted  to  address  as  one  of  our  experiments  during  the  campaign  was 
to  collect  data  from  two  different  receiving  geometries  during  the  same  satellite  pass.  The  goals 
of  this  are  to  (1)  have  at  least  one  observation  not  viewing  the  heated  region  with  a  geometry 
parallel  to  the  local  geomagnetic  field  (avoid  the  propagation  enhancement  to  scintillation 
observed  when  viewing  along  the  local  field  direction),  and  (2)  sample  different  regions  of  the 
heated  region  to  look  at  differences  in  the  effects  at,  and  away  from,  the  magnetic  zenith  during 
the  same  heating  experiment.  Our  experiment  on  24  February  2008  using  the  Oscar  25  satellite 
was  designed  to  look  at  these  issues.  We  deployed  our  Mobile  receiver  away  from  the  HAARP 
site  to  provide  a  different  viewing  geometry  than  that  available  from  the  Gakona  ITS  receiver, 
and  we  used  a  new  beam  geometry  developed  by  Dr.  K.  Groves  of  AFRL  and  Dr.  M.  McCarrick 
of  BAE  Systems.  This  new  geometry  spreads  the  beam  in  a  fan  shape  with  the  beam  center 
aimed  at  the  local  zenith  (not  magnetic)  and  the  fan  semi-major  axis  aligned  with  the  array 
orientation  (just  east  of  geographic  north). 

Figure  27  shows  the  geometry  for  this  pass  and  the  resulting  VHF  intensity  and  VHF/UHF 
differential  phase  data  collected  during  this  pass.  From  the  IPP  track  map  in  Figure  27a,  you  can 
see  that  the  IPP  from  the  Mobile  receiver  (shown  in  blue)  passed  through  the  southern  end  of  the 
heated  region  and  sampled  very  near  the  IRI  magnetic  zenith  position  which  the  Gakona  receiver 
geometry  (shown  in  red)  samples  the  northern  end  of  the  heated  region  well  away  from  the  IRI 
magnetic  zenith.  The  elongated  ovoids  shown  in  the  map  is  the  mapping  of  the  IRI  beam  3-dB 
and  6-dB  contours  to  a  220-km  altitude.  The  times  that  the  ray  paths  were  within  these  contours 
are  indicated  in  the  data  plots  in  Figure  27b  by  horizontal  heavy  (3-dB)  and  light  (6-dB)  lines. 
The  heavy  vertical  line  in  these  plots  indicates  where  the  geometry  crosses  the  semi-major  axis 
of  the  heated  region.  Note  that  this  pass  proceeded  from  south  to  north. 

The  results  seen  in  the  data  from  the  Mobile  pass  (top  two  panels  in  Figure  27b)  are  similar 
to  those  observed  in  previous  IRI  experiments  using  signals  from  these  LEO  satellite  beacons. 
The  results  from  the  Gakona  receiver  are  different  and  quite  interesting.  In  particular,  the  most 
disturbed  section  of  the  pass  occurs  not  when  the  signal  was  passing  through  the  heated  region, 
but  to  the  south  of  that  region  at  a  similar  latitude  to  where  the  effects  were  seen  in  the  data  from 
the  Mobile  receiver.  We  have  not  fully  analyzed  this  pass  yet,  but  one  initial  interpretation  is 
that  what  we  are  seeing  is  structure  that  has  been  advected  out  of  the  heated  region  by 
background  flow  of  the  ionosphere  (which  is  supported  by  preliminary  analyses  of  the  Digisonde 
drift  data)  with  a  possible  suggestion  from  the  magnitude  of  the  enhancement  seen  in  the  Gakona 
phase  data  that  these  structures  may  be  elongated  across  magnetic  field  lines  and  along  the  local 
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L-shell,  and  not  just  along  the  magnetic  field.  This  is  a  very  preliminary  finding  which  we  plan 
to  pursue  during  the  next  few  months. 
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Figure  27.  (a)  220-km  Ionospheric  Penetration  Point  (IPP)  track  for  the  24  Feb  2008  Oscar  23  pass  for  the 
raypath  from  the  Gakona  receiver  (in  red)  and  the  NWRA  Mobile  receiver  (in  blue).  The  extended 
ovoid  shapes  indicate  the  3-dB  and  6-dB  contours  of  the  IRI  beam  used  in  this  experiment,  (b) 
Detrended  (30-second  cutoff)  VHF  intensity  and  phase  as  observed  from  the  NWRA  Mobile  receiver 
(top  two  panels),  and  as  observed  from  the  Gakona  receiver  (bottom  two  panels.  Heavy  vertical  line 
indicates  closest  approach  to  the  IRI  beam  semi-major  axis,  heavy  horizontal  line  indicates  where  the 
IPP  was  within  the  3-dB  contours,  and  the  light  horizontal  line  indicates  where  the  IPP  was  within  the  6- 
dB  contours. 

2.5.23.  GPS  Observations 

Dr.  Gennady  Milikh,  University  of  Maryland,  ran  a  series  of  experiments  with  the  HAARP 
IRI  that  employed  the  NWRA  GPS  receiver  located  at  the  HAARP  site  as  a  primary  diagnostic. 
Figure  28  shows  the  Total  Electron  Content  (TEC)  derived  from  one  of  Dr.  Milikh’s  experiment 
periods  showing  clear  evidence  of  structures  generated  by  the  HAARP  IRI  during  the 
experiment.  We  have  submitted  a  paper  to  Geophysical  Research  Letters  describing  this  result 
[Milikh  et  al.,  2008]. 
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Figure  28.  Data  from  the  NWRA  HAARP  GPS  receiver  taken  during  an  experiment  by  Dr.  G.  Milikh.  The 
bottom  panel  is  the  absolute  slant-TEC  derived  from  observations  of  the  time-delay  and  phase-advance 
on  the  signals  from  GPS  satellite  PRN  21,  and  the  top  panel  is  the  data  detrended  by  a  high-pass  filter 
with  a  40  second  cut-off  period.  The  horizontal  dashed  line  indicated  times  that  the  IRI  transmitter  was 
on  and  off  (10-second  intervals). 

3.  NON-HAARP  TOPICS  -  TECHNICAL  AREA  3(A) 

Contract  work  on  topics  under  BAA  Technical  Area  3(a)  during  the  period  covered  by  this 
report  was  focused  on  two  primary  areas  of  interest:  (1)  ionospheric  scintillation,  and  (2) 
ionospheric  Total  Electron  Content  (TEC)  derived  from  measurements  using  GPS  signals. 

3.1.  Ionospheric  Scintillation  Research 

Two  research  efforts  were  undertaken  during  the  contract  period  focused  on  ionospheric 
scintillation.  Continuing  the  thrust  of  the  work  reported  in  SRI  (Section  3.1.)  wherein 
algorithms  for  specifying  the  uncertainty  in  estimates  of  scintillation  indices  calculated  by  the 
WBMOD  ionospheric  scintillation  model  [ Secan ,  2004]  were  developed  and  added  to  that 
model,  we  used  several  years  of  scintillation  observations  collected  by  the  ITS  1  OS  receiver 
located  at  the  HAARP  site  in  Gakona,  AK,  to  assess  the  performance  of  the  WBMOD  model  at 
an  auroral  station.  Results  of  that  work  were  provided  in  a  separate  report  (Secan,  2005). 

The  second  effort  in  this  area  consisted  of  two  studies  looking  at  ionospheric  scintillation 
models  and  assessed  their  performance  in  terms  of  scintillation  observations  made  by  several 
systems  which  are  operationally  impacted  by  the  effects  of  scintillation.  The  results  of  these 
studies  are  presented  in  Appendix  D  to  this  report. 


43 


3.2.  Use  of  GPS  Signals  to  Derive  Ionospheric  TEC 

The  use  of  radio  beacons  on  satellites  to  measure  the  Total  Electron  Content  (TEC)  between 
the  ground  and  the  satellite  dates  back  to  the  first  artificial  satellite,  Sputnik  I,  in  1957  {Al’pert  et 
al.,  1958).  This  parameter  became  of  interest  to  the  DoD  community  as  TEC  provides  a  useful 
observation  for  estimating  the  impact  of  the  ionosphere  on  a  wide  range  of  communications  and 
surveillance  systems  which  make  use  of  the  transionospheric  propagation  channel.  In  the  1970s 
and  well  into  the  1990s,  a  primary  source  of  ionospheric  TEC  for  operational  use  was  a  network 
of  polarimeters,  which  make  use  of  the  Faraday  rotation  effect,  established  by  the  U.S.  Air  Force 
Cambridge  Research  Laboratory  (AFCRL)  for  use  in  models  run  by  the  USAF  Air  Weather 
Service  (AWS)  ( Mendillo  and  Klobuchar,  1974;  Flattery  and  Ramsay,  1975).  With  the  launch 
of  the  first  operational  Global  Positioning  System  (GPS)  satellites  in  the  late  1980s,  work  was 
already  underway  to  use  the  two  phase-coherent  GPS  L-band  signals  (LI  and  L2)  as  a 
replacement  for  the  aging  polarimeter  network  (, Klobuchar  and  Rino,  1983).  This  work 
culminated  in  the  mid  1990s  with  development  and  deployment  of  the  USAF  Ionospheric 
Measuring  System  (IMS)  GPS-based  TEC  systems  {Bishop  et  al.,  1994). 

While  GPS-based  TEC  was  a  major  advance  over  the  old  polarimeter  network  TEC,  there 
remained  two  issues  to  resolve;  one  that  was  a  carry  over  of  an  unresolved  issue  with  the 
polarimeter-based  TEC  and  a  second  due  to  limitations  of  the  polarimeter-based  TEC.  The  first 
issue  is  how  to  use  TEC,  an  integral  measure,  to  updated  models  of  the  full  electron  density 
profile.  This  will  be  an  issue  for  any  system  that  measures  TEC.  The  second  issue  is  a 
somewhat  artificial  one,  although  important  in  an  operational  setting.  The  polarimeter-based 
TEC  was  able  to  measure  TEC  only  to  an  altitude  of  1,500  to  2,000  km  due  to  the  weighting 
effect  of  the  earth’s  magnetic  field  on  the  Faraday  rotation  phenomenon  on  which  these  TEC 
measurements  were  made  {Titheridge,  1972).  This  limitation  led  to  an  artificial  “cap”  on  the 
upper  altitude  range  of  models  used  by  the  USAF  AWS;  a  limitation  that  carried  forward  to 
newer  models  and  that  has  led  to  a  need  for  estimates  of  TEC  up  to  an  altitude  of  2,000  km.  This 
has  been  codified  into  a  specification  of  “ionospheric  TEC”  to  refer  to  TEC  up  to  this  altitude 
and  “plasmaspheric  TEC”  for  TEC  from  that  altitude  out  to  infinity. 

The  issue  of  the  impact  of  plasmaspheric  TEC  also  has  ramifications  on  the  ability  to 
convert  TEC  observed  to  GPS  satellites  from  a  measure  along  a  raypath  of  any  orientation 
(called  slant  TEC)  to  an  estimate  of  the  equivalent- vertical  TEC.  This  was  also  a  parameter  used 
by  older  AWS  models,  and  the  process  of  converting  from  slant  TEC  to  equivalent- vertical  TEC 
is  a  key  part  of  the  process  developed  for  the  IMS  system  to  estimate  the  biases  in  differential 
group  delay  (DGD)  estimates  of  TEC  due  to  uncalibrated  delays  (between  the  LI  and  L2  signals) 
in  both  satellites  and  receivers  {Bishop  et  al.,  1995),  known  as  the  SCORE  process  (for  Self 
Calibration  Of  Range  Errors).  The  effect  of  plasmaspheric  TEC  on  this  process  is  described  in 
Lunt  et  al.,  [1999],  A  preliminary  method  for  dealing  with  this  problem,  which  also  addressed 
the  issue  of  separating  the  TEC  into  ionospheric  and  plasmaspheric  components,  was  developed 
by  NWRA  and  denoted  SCORPION  (SCOre  for  Plasmasphere  and  IONosphere)  {Mazzella  et  al., 
2002). 

The  original  plan  for  this  effort  was  to  continue  development  of  the  SCORPION  algorithm 
as  a  way  to  resolve  at  least  the  issue  of  partitioning  GPS  slant  TEC  into  an  equivalent-vertical 
ionospheric  TEC  and  a  plasmaspheric  TEC.  That  effort  was  funded  by  AFRL  via  a  Small 
Business  Innovative  Research  contract  {Rao,  2004;  Rao  and  Mazzella,  2007).  To  complement 
that  separate  effort,  and  to  address  issues  other  than  the  partitioning  of  TEC  into  ionospheric  and 
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plasmaspheric  components,  three  parallel  efforts  were  conducted  under  this  contract  looking  at 
improving  use  of  GPS  TEC  data: 

1.  Expanding  on  an  existing  data-assimilation  system,  the  GPS  Ionospheric  Inversion  (GPSII) 
program  ( Fridman  and  Nickisch,  2001),  to  deal  with  the  issue  of  converting  GPS  slant  TEC 
measurements  into  electron  density  profdes  and,  as  a  side  effect,  with  the  partitioning  issue  as 
well. 

2.  Expand  the  existing  SCORE  processing  to  look  at  data  from  a  chain  of  stations  along  a 
longitude  sector  to  assess  the  performance  of  SCORE  in  dealing  with  the  impact  of  the 
plasmaspheric  TEC  during  solar  minimum  conditions. 

3.  Begin  work  to  codify  and  quantify  the  various  sources  of  uncertainty  in  GPS  TEC 
measurements  (both  slant  and  equivalent-vertical). 

4.  Deploy  GPS  systems  to  a  high  latitude  station  (Clear  AFS,  AK)  to  run  the  SCORPION 
analysis  software  to  both  test  and  demonstrate  that  the  SCORPION  system  provided  realistic  bias 
results  at  a  location  where  there  should  be  little,  or  no,  contribution  from  the  plasmasphere. 

We  were  fortunate  in  that  there  were  a  number  of  opportunities  to  test  and  exercise  the  new 
capabilities  developed  under  these  efforts  to  meet  the  needs  of  AF  users.  Appendix  E  describes 
extensions  made  to  GPSII  to  better  use  GPS  TEC  data  and  several  applications  of  this  new 
technology  to  specific  applications.  We  were  also  able  to  make  use  of  the  GPSII  developments 
to  provide  an  extensive  analysis  of  data  in  the  west  European  longitude  sector  in  support  of 
analysis  by  a  DoD  agency  of  the  performance  of  radar  systems  in  the  UK  (the  spacetrack  radar  at 
Fylingdales).  The  work  on  assessing  the  performance  of  SCORE  from  a  chain  of  stations  along 
the  west  cost  of  the  U.S.  and  our  work  to  quantify  uncertainty  in  GPS  TEC  measurements  were 
leveraged  to  provide  support  for  a  validation  and  calibration  effort  for  the  DMSP  SSUSI 
ultraviolet  spectrometer  sensor.  In  this  section,  we  will  provide  an  overview  of  the  SCORE 
system  developed  for  our  research,  describe  the  two  efforts  just  mentioned,  and  the  data- 
collection  deployment  to  Clear  AFS,  in  the  context  of  the  core  focus  of  the  overall  effort  to 
improve  the  use  of  GPS  TEC  measurements. 

3.2.1.  SCORE  Test  and  Analysis  Implementation 

As  a  precursor  to  the  research  planned  that  would  make  use  of  the  SCORE  analysis  system, 
we  spent  some  time  developing  a  more  flexible  implementation  of  the  SCORE  process  that 
would  provide  a  more  user-friendly  and  extensible  platform  for  use  not  only  in  our  planned 
research,  but  also  for  future  use.  The  existing  SCORE  implementation  was  tied  strongly  to  the 
real-time  implementation  used  in  the  IMS  system,  and  the  test  implementation  that  had  been 
derived  from  the  operational  implementation  was  complicated  to  set  up  and  cumbersome  to  use. 
There  was  also  very  limited  information  available  upon  which  to  base  assessments  of  how  well 
the  analysis  worked  and  the  uncertainty  in  the  TEC  generated  by  the  analysis.  The  existing 
implementation  was  also  tied  to  the  Windows  operating  system,  and  we  wanted  a  system  that 
could  run  on  a  wide  range  of  platforms. 

The  decision  was  made  to  build  the  analysis  system  around  the  existing  Fortran  codes  that 
make  up  the  SCORE  analysis  system,  making  only  those  changes  to  these  codes  as  necessary  to 
add  capabilities  not  currently  in  existence.  The  scripting  codes  that  were  built  to  run  these 
programs,  however,  were  to  be  completely  rewritten.  For  the  sake  of  portability,  we  decided  to 
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use  the  python  scripting  language  for  the  majority  of  the  processing  scripts  (wherein  lies  most  of 
the  expertise  in  running  the  SCORE  Fortran  programs)  and  the  Generic  Mapping  Tools  (GMT) 
package  for  all  display  generation.  Both  of  these  are  freely  available  over  the  Internet,  and 
provide  support  on  all  major  computing  platforms.  The  system  was  also  designed  to  present  a 
range  of  graphical  outputs  in  a  form  that  is  easily  accessible  by  standard  Internet  browsers  such 
as  Microsoft’s  Internet  Explorer  and  the  Mozilla  Firefox  browser.  One  final  constraint  on  the 
system  was  that  it  would  be  designed  such  that  replacing  the  SCORE  bias-determination  process 
by  the  SCORPION  process  should  be  relatively  simple.  For  our  needs,  we  tested  the  new  codes 
on  Windows  (XP),  Linux,  and  Unix  (Solaris)  systems.  As  a  starting  point  for  the  core  processing 
script,  we  used  a  script  that  had  been  developed  by  S.  Rao  for  testing  the  SCORPION  system 
[Rao,  2004], 

Details  for  this  rather  elaborate  system,  which  covers  the  processing  from  downloading 
RINEX  data  files  from  sites  around  the  Internet  to  display  of  both  TEC  data  and  a  range  of 
diagnostics  and  tools  for  repackaging  the  TEC  data  for  delivery  to  other  users,  are  described  in 
detail  in  formal  documentation  written  for  this  system  and  delivered  to  AFRL  ( Secan  and 
Holland,  2008).  This  report  will  describe  this  system  at  a  high  functional  level,  focusing  mostly 
on  the  new  tools  developed  for  assessing  performance  and  in  ways  it  has  been  used  in  two 
capabilities  demonstrations.  Examples  of  the  primary  displays  that  this  system  is  capable  of 
producing  are  presented  in  Appendix  F  to  this  report. 

The  system  is  designed  such  that  it  can  be  used  to  download  RINEX- format  data  from  any 
repository  of  RINEX  data  files  on  the  Internet,  process  the  data  using  the  SCORE  system, 
generate  a  number  of  displays  from  the  resulting  analysis,  store  the  results  of  the  analysis  in  an 
archive  database  system,  and  post  the  displays  in  a  cross-linked  web  site  for  ease  in  viewing. 
This  process  can  be  configured  to  do  one  day  for  a  single  station  or  multiple  days  for  a  collection 
of  stations.  A  number  of  tailored  versions  of  this  system  were  developed  in  the  course  of  the 
research  conducted  on  this  contract,  ranging  from  a  hands-on  one  station-day  per  run  setup,  to  a 
fully  automated  system  that  downloads,  processes,  and  stores  data  from  a  chain  of  stations  on  a 
daily  basis  without  human  intervention. 

The  tools  developed  to  aid  in  assessing  the  performance  of  the  bias-calibration  process  and 
the  uncertainties  in  the  TEC  derived  based  on  those  biases  can  be  categorized  as  follows: 

1 .  Displays  of  the  resulting  equivalent- vertical  TEC  that  are  simple  to  generate  and  easy  to  use. 

2.  Statistics  generated  from  the  TEC  produced  by  the  process 

3.  Comparison  of  the  biases  generated  by  the  SCORE  (or  SCORPION)  process  with  biases  from 
other  days  or  from  other  sources. 

4.  Comparison  of  the  TEC  generated  with  TEC  from  other  sources. 

Looking  at  the  display  examples  shown  in  Appendix  F,  Figures  FI  and  F2  are  in  the  first 
category;  parts  of  Figures  FI  and  F2  (the  average  and  standard  deviation  parameters  included  in 
these  plots)  and  Figure  F4  are  in  the  second  category;  Figures  F3,  F8,  and  F9  are  in  the  third 
category,  and  Figures  F5,  F6,  and  F7  are  in  the  fourth  category.  The  set  of  Figures  FI  through 
F4  can  be  generated  routinely  and  posted  to  a  cross-linked  web  site  for  review,  as  can  Figures  F5 
through  F6  if  data  are  available  from  the  JASON  spacecraft  (see  Section  3.2.3.).  The  other 
figures  can  be  generated  automatically,  but  they  are  set  up  to  be  used  by  what  we  refer  to  as  a 
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“skilled  analyst”  to  use  in  assessing  how  well  or  poorly  SCORE  (or  SCORPION)  has  performed 
on  a  particular  data  set. 

This  concept  of  the  “skilled  analyst”  is  a  stepping  stone  to  a  more  automated  system  for 
assessing  overall  performance.  The  plan  is  that  by  using  these  tools  manually  to  determine  how 
the  system  is  behaving,  we  will  be  able  to  identify  which  of  the  various  metrics  we  have 
established  for  assessing  system  performance  could  be  worked  into  an  automated  system.  At  this 
time,  the  expertise  in  this  area  is  at  the  level  of  a  person  who  understands  both  the  SCORE 
system  and  the  behaviors  and  limitations  of  GPS  TEC  data,  who  can  look  at  the  various  displays 
and  output  parameters  and  make  a  judgment  as  to  whether  the  biases  generated  for  a  given  day, 
and  the  resulting  TEC,  are  good  or  not.  Some  of  the  outputs,  such  as  the  average  TEC  as  a 
function  of  local  solar  time  and  the  standard  deviation  of  the  individual  TEC  observations  about 
this  mean,  can  be  used  as  preliminary  assessments  of  quality  and  uncertainty,  and  have  been 
provided  as  such  to  end  users  (see  Section  3.2.3.).  Further  automation  of  this  process  will 
require  further  work  with  these  tools. 

3.2.2.  Western  Europe  Sector  GPSII  Analysis 

An  earlier  demonstration  test  of  an  early  version  of  the  SCORE  implementation  just 
described  and  of  the  GPSII  system  described  in  Appendix  E  was  afforded  by  a  request  from  a 
USAF  agency  for  analysis  of  data  from  several  time  periods  around  solar  maximum  to  be  used  in 
assessing  the  performance  of  new  software  developed  for  early  warning  radars  to  mitigate  the 
impact  of  the  ionosphere  on  radar  operations.  By  the  end  of  this  effort,  we  provided  estimates  of 
ionospheric  TEC  along  user-specified  ground-to-satellite  ray  paths  covering  three  two-week 
periods:  05  to  19  January  2002,  09  to  23  June  2001,  and  27  March  to  10  April  2001.  For  each 
day  of  these  two-week  periods,  RINEX  data  for  a  set  of  stations  in  a  region  centered  over  the 
United  Kingdom  and  reaching  into  northern  Africa  on  the  south,  eastern  Europe  on  the  east, 
Scandinavia  on  the  north,  and  southern  Greenland  on  the  west,  was  downloaded  and  processed 
by  the  SCORE  implementation  to  provide  an  initial  analysis  and  then,  by  the  GPSII  system,  to 
provide  grids  of  electron  density  profiles  from  which  TEC  along  any  trajectory  within  the  grid 
could  be  easily  calculated. 

While  the  SCORE  output  was  not  used  directly  in  generating  the  end  product  (TEC  along 
specific  trajectories),  it  was  used  to  identify  stations  that  had  data-quality  issues  and  as  a  cross¬ 
check  to  the  TEC  produced  by  the  GPSII  analysis  to  insure  that  the  newer  technology  in  the 
GPSII  software  was  producing  estimates  of  satellite  and  receiver  biases  that  were  consistent  with 
those  estimated  by  the  SCORE  process.  We  used  SCORE  as  a  pre-selector  tool  because  it  ran 
much  quicker  than  GPSII  and  works  with  one  station  independent  of  all  other  data.  This  allowed 
us  to  make  a  quick-look  assessment  of  the  data  to  be  used  with  no  possibility  that  odd  behavior  at 
any  one  given  station  was  being  caused  by  data  problems  at  some  other  station  that  was  part  of 
the  same  analysis.  The  cross-check  of  the  biases  generated  by  the  two  approaches  provided  a 
“sanity  check”  on  the  more  complex  GPSII  analysis. 

Plots  of  the  daily  TEC  similar  to  those  shown  in  Figures  FI  and  F2  were  generated  for  both 
the  SCORE  analysis  and,  when  the  analysis  output  became  available,  for  the  GPSII  analysis. 
The  GPSII  software  was  modified  to  provide,  as  a  diagnostic  output,  the  TEC  produced  by  the 
“forward  model”  within  GPSII  for  GPS  TEC  corresponding  to  each  input  GPS  TEC.  If  the  two 
processes  were  generating  completely  consistent  estimates  for  TEC,  the  plots  would  be  identical. 
Two  differences  in  how  these  systems  operate  make  an  identical  analysis  unlikely.  First,  the 
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SCORE  process  converts  the  observed  slant  TEC  into  an  equivalent-vertical  TEC,  while  the 
GPSII  process  generates  the  vertical  TEC  estimates  by  a  true  vertical  integration,  thus,  there  are 
differences  due  to  known  deficiencies  in  the  process  used  within  SCORE.  Second,  the  SCORE 
process  analyzes  the  data  one  station  at  a  time,  with  no  attempt  made  to  reconcile  the  results  with 
results  from  other  stations,  while  GPSII  attempts  to  find  solutions  that  provide  the  best  possible 
fit  to  all  stations  for  which  it  has  data.  The  SCORE  plots  were  used  initially  in  the  input-data  QC 
process  in  which  we  eliminated  unsuitable  stations.  After  the  GPSII  data  were  available,  cross¬ 
comparisons  were  made  to  insure  that  the  GPSII  analysis  was  generating  realistic  grids  of 
ionospheric  electron  density  in  the  analysis  region. 

All  trajectory  data  were  generated  and  provided  to  the  USAF  customer  by  the  end  of  May 
2007.  This  exercise  assisted  us  in  refining  several  of  the  tools  we  had  developed  for  the  SCORE 
implementation,  and  led  to  several  improvements  to  the  GPSII  software,  as  described  in 
Appendix  E.  One  major  improvement  was  in  the  runtime  of  the  GPSII  code,  which  was 
decreased  by  over  a  factor  of  two  in  order  to  handle  the  large  volume  of  data  employed  in  this 
analysis. 

3.2.3.  SCORPION  Performance  at  Auroral  Latitudes 

Shortly  after  the  study  described  above,  we  were  provided  an  opportunity  to  deploy  a  GPS 
receiver  system  operating  the  NWRA  GIOS  (GPS  Ionospheric  Observing  System)  software 
package  to  Clear  AFS,  AK,  in  support  of  a  USAF  test  activity,  Combat  SPACETRACK 
Integration  Node  Global  Enhanced  Radar  (STINGER).  Rather  than  the  usual  SCORE  calibration 
package,  the  system  was  deployed  with  a  real-time  version  of  the  SCORPION  calibration 
package  which  had  just  been  completed  under  an  SBIR  contract  ( Rao  and  Mazzella,  2007).  The 
system  was  deployed  to  Clear  AFS  in  mid- July  2007,  and  installed  and  operated  by  NWRA 
consultant.  Mr.  Frank  Smith,  of  Northwest  Research  Engineering  (NRE).  Data  was  collected  in 
conjunction  with  tests  run  using  the  radar  system  located  at  Clear  AFS  as  part  of  the  Combat 
STINGER  testing  for  the  period  24  through  29  July  2007. 

The  GPS  TEC  data  collected  was  calibrated  in  real-time  using  the  SCORPION  software  run 
in  a  mode  designed  to  replicate  the  results  that  would  be  obtained  from  SCORE.  This  run  mode 
was  designed  to  be  used  when  the  data  were  being  collected  at  a  location  far  enough  poleward  of 
the  plasmasphere  boundary  that  there  should  be  no,  or  minimal,  impact  by  the  plasmasphere  on 
the  calibration  operation.  This  deployment  provided  an  opportunity  to  collect  data  at  such  a 
location  to  test  the  real-time  operation  of  SCORPION  using  this  run  mode. 

The  performance  of  SCORPION  was  assessed  by  comparing  the  estimates  of  the  combined 
receiver  and  satellite  biases  obtained  by  SCORPION  in  the  field  with  estimates  generated  by  the 
production  version  of  SCORE  analyzing  the  data  after-the-fact.  The  processing  system 
developed  for  testing  SCORE  performance,  as  described  in  the  previous  section,  modified  to  use 
data  from  the  GIOS  system  rather  than  RINEX-format  data,  was  used  for  this  testing.  For  the 
six-day  period  of  the  test,  the  RMS  difference  between  the  SCORPION  and  SCORE  estimates  of 
the  biases  was  0.10  TECU,  and  all  were  within  1  TECU.  The  majority  of  the  large  differences 
occurred  on  the  initial  day  of  the  data  collection,  and  appear  to  have  been  due  to  problems 
getting  the  data  collection  process  started.  These  problems  were  not  reflected  in  the  SCORE 
processing,  so  the  two  analysis  systems  were  working  with  slightly  different  data  sets  for  that 
day.  Overall,  the  comparison  showed  that  SCORPION  operating  in  a  SCORE  mode  (ionosphere 
only)  performs  as  well  as  the  production  version  of  SCORE.  This  was  the  desired  outcome,  as 


48 


SCORE  has  proven  to  be  as  good  a  process  as  any  available  for  determining  the  biases  from  a 
single  station  in  near  real-time. 

The  Combat  STINGER  program  needed  to  have  these  observations  in  a  form  that  was  more 
easily  used  in  their  post  analysis  software.  We  obtained  the  services  of  a  consultant,  Dr.  Rob 
Daniell,  to  take  the  data  we  had  collected  and  provide  Combat  STINGER  with  the  analyses  in 
the  format  they  needed.  These  data  were  delivered  to  AFRL  for  delivery  to  the  Combat 
STINGER  team  in  September  2007. 

3.2.4.  SCORE  West-Coast  US  Analysis 

The  SCORE  implementation  reached  its  full  maturity  in  a  final  application  in  which  RINEX 
data  from  a  chain  of  stations  along  the  west  coast  of  the  United  States  was  analyzed  using 
SCORE  to  provide  estimates  of  TEC  for  use  in  a  calibration  and  validation  (cal/val)  study  of 
TEC  data  produced  by  the  new  DMSP  Special  Sensor  Ultraviolet  Spectrographic  Imager 
(SSUSI).  We  used  this  opportunity  to  extend  the  tools  already  developed  and  to  develop 
additional  tools  to  allow  us  to  compare  TEC  derived  from  the  SCORE  system  with  TEC  from 
another  source,  the  JASON  oceanic  altimeter  data.  Many  of  the  graphical  tools  shown  in 
Appendix  F  were  developed  for  this  effort,  which  has  greatly  expanded  the  toolkit  available  for 
assessing  both  the  performance  of  the  SCORE  system  and  quantifying  the  uncertainty  in  TEC 
derived  by  that  process.  Table  4  lists  the  nine  stations  for  which  data  were  provided,  and  Table  5 
lists  the  days  in  2007  for  which  slant-range  TEC  data  were  provided  for  each  of  the  nine  stations 
in  Table  4.  All  data  identified  by  these  two  tables  were  provided  to  AFRL  and  Aerospace  during 
the  report  period. 


Table  4.  List  of  stations  for  which  slant-range  TEC  data  derived  from  GPS  observations  were  provided. 


Station 

Geographic 

Geomagnetic 

Mlat  Difference 

Name 

ID 

Lat 

Lon 

Lat 

Lon 

To  South 

To  North 

Neah  Bay 

neah 

48.30 

235.38 

54.08 

293.93 

1.94 

N/A 

Fort  Stevens 

ftsl 

46.21 

236.04 

52.14 

295.43 

1.48 

1.94 

Corvalis 

corv 

44.59 

236.69 

50.66 

296.70 

2.75 

1.48 

Yreka 

ybhb 

41.73 

237.29 

47.91 

298.28 

3.49 

2.75 

Petaluma 

pl98 

38.26 

237.39 

44.42 

299.47 

2.48 

3.49 

Bakersfield 

bkrl 

35.13 

240.89 

41.94 

304.09 

2.17 

2.48 

Point  Loma 

plo5 

32.67 

242.76 

39.77 

306.72 

7.61 

2.17 

La  Paz 

lpaz 

24.14 

249.68 

32.16 

315.89 

4.58 

7.61 

Colima 

col2 

19.24 

256.30 

27.58 

323.40 

N/A 

4.58 

Table  5.  Dates  for  which  slant-TEC  data  derived  from  the  chain  of  stations  shown  in  Table  4  were  provided. 


#  Days 

Data  Days  (2007) 

Start 

End 

40 

30-Mar-07 

08-May-07 

11 

21 -May-07 

31 -May-07 

14 

12-Jul-07 

25-Jul-07 

21 

10-Sep-07 

30-Sep-07 

12 

15-Dec-07 

26-Dec-07 
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Figures  F5  through  F7  show  a  set  of  analysis  displays  that  were  developed  to  facilitate  the 
comparison  of  GPS  TEC  derived  using  SCORE  and  TEC  derived  from  the  JASON  data.  Figure 
F5  is  a  composite  plot,  so  named  because  the  GPS  TEC  shown  is  a  composite  of  GPS  TEC  data 
averaged  to  provide  a  comparison  with  the  JASON  data.  This  composite  average  (the  line  of 
black  error  bars)  is  generated  by  selecting  those  GPS  TEC  values  that  are  within  a  user-selected 
spatial  and  temporal  distance  from  the  track  along  which  JASON  data  are  available  (the  line  of 
solid  blue  dots)  and  creating  a  weighted  average,  where  the  weighting  is  dependent  on  the 
distance  (spatial  and  temporal)  between  the  GPS  TEC  observation  and  the  JASON  TEC 
observation.  The  width  of  the  error  bars  on  the  GPS  TEC  points  is  given  by  the  weighted 
standard  deviation  of  the  GPS  TEC  data  used  to  generate  the  average.  The  green  points  along 
the  bottom  of  this  plot  indicate  the  difference  between  the  GPS  composite  TEC  and  the  JASON 
TEC.  Note  that  previous  analyses  of  TEC  estimates  calculated  from  the  JASON  altimeter  data 
found  that  JASON-derived  TEC  overestimates  the  actual  TEC  by  4-5  TECU  ( Bishop ,  private 
communication) . 

Figures  F6  and  F7  provide  an  in-depth  look  at  the  individual  GPS  TEC  data  that  were  used 
in  generating  the  average  GPS  TEC  shown  on  Figure  F6.  The  difference  between  these  two 
figures  is  the  plot  in  the  lower-right  quadrant  of  the  figure.  In  Figure  F6,  the  location  of  the  GPS 
TEC  values  and  the  location  of  the  JASON  track  are  shown  as  a  function  of  azimuth  and 
elevation  angle  from  the  site,  in  this  case  St.  John’s,  Newfoundland.  In  Figure  F7,  this  is 
replaced  by  a  plot  of  the  locations  in  latitude  and  longitude,  with  the  sub-satellite  track  location 
used  for  the  JASON  data  and  the  350-km  Ionospheric  Penetration  Point  (IPP)  used  for  the  GPS 
TEC  data.  The  purpose  of  these  plots  is  to  provide  an  analyst  a  more  detailed  look  at  the 
locations  of  the  GPS  TEC  data  with  respect  to  the  JASON  TEC  locations  (in  space  and  time)  and 
of  the  scatter  of  the  GPS  TEC  as  a  function  of  local  time  and  latitude. 

A  detailed  look  at  some  of  the  results  derived  from  this  effort,  including  a  multi-station 
version  of  the  plot  shown  in  Figure  F5,  is  presented  in  Bishop  et  al.  [2008],  This  preliminary 
look  at  the  performance  of  SCORE  and  the  impact  of  the  plasmasphere  during  solar  minimum 
conditions  showed  that  these  new  tools  will  be  useful  for  assessing  the  overall  performance  of 
the  SCORE  process  and,  eventually,  the  SCORPION  process. 
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Appendix  A 

Tomography  Image  Coverage 

Tables  A1  through  A5  provide  a  summary  of  images  generated  from  the  Alaskan 
tomography  network  from  2004  through  2008  (31  July),  respectively,  as  a  function  of  month. 
The  upper  panel  separates  the  images  by  the  set  of  stations  used  to  generate  each  image.  The 
entries  in  the  columns  under  "Stations"  indicate  which  stations  contributed  to  the  image;  for 
example,  the  top  row  lists  images  generated  from  the  station-set  Kaktovik,  Fort  Yukon,  and 
Delta  Junction.  The  bottom  panel  shows  the  percentage  of  images  in  each  month  that  included 
data  from  the  station  listed  in  the  left-most  column. 

Note  that  the  decrease  in  the  number  of  images  and  in  the  number  of  stations  from  the  chain 
in  these  images  towards  the  end  of  the  report  period  were  due  to  two  factors:  (1)  the  primary 
satellites  used  for  the  tomography  work,  the  Navy  Ionospheric  Measuring  System  (NIMS),  or 
Oscar,  satellites,  are  beginning  to  go  off  the  air;  and  (b)  the  two  UAF  GI  stations,  Kaktovik  and 
Delta  Junction,  became  highly  unreliable  and  then  left  the  air  completely  during  2007  and  2008. 
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Table  A2.  Station-Availability  Summary  for  Tomographic  Images  in  2005 
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Table  A3.  Station- Availability  Summary  for  Tomographic  Images  in  2006 
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Table  A4.  Station-Availability  Summary  for  Tomographic  Images  in  2007 
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Table  A5.  Station-Availability  Summary  for  Tomographic  Images  in  2008  (Through  July) 


Number  of  Images 
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Appendix  B 

ITS  Receiver  Stations  Status 

The  following  summary  report  of  the  status  of  the  NWRA-operated  ITS  tomography 
receiver  chain  was  provided  by  NWRA  consultant  Mr.  Frank  Smith.  NWRA  maintains  and 
operates  five  ITS  1  OS  Receiver  systems  for  the  HAARP  tomography  chain.  At  any  time,  four  of 
these  systems  are  deployed  at  sites  in  Alaska,  while  the  fifth  system  is  a  reserve  to  serve  as  back¬ 
up  and  a  test  platform  for  software  or  hardware  upgrades  or  trouble-shooting. 

Bl.  Description  of  ITS  Receiver  Systems 

The  five  ITS  10  receiver  systems  are  nearly  identical,  comprising  a  two-band  coherent 
receiver  (VHF  and  UHF),  an  associated  two-band  antenna  and  a  computer.  The  computer 
controls  the  operation  of  and  collects  data  from  the  receiver,  and  it,  in  turn,  is  remotely  accesses 
to  install  new  software,  update  the  NORAD  Two  Line  Elements  (TLE)  and  adjust  various 
operational  parameters. 

The  VHF  Detector  section  of  the  receiver  is  comprised  almost  entirely  of  obsolete  parts. 
However,  over  the  years,  we  have  had  few  failures  with  this  section,  but  the  replacement  parts 
are  nearly  used  up.  Another  part  of  the  receiver,  the  Synthesizer,  is  also  made  with  an  obsolete 
oscillator  for  which  no  replacement  is  available.  These  oscillators,  however,  have  never  failed. 
The  remainder  of  the  receiver,  while  comprised  of  1990's  era  parts,  remains  serviceable. 

The  computers  are  quite  old  and  obsolete:  300  MHz  machines  running  Windows  NT. 
Additionally,  of  the  two  National  Instruments  cards  that  interface  between  the  receiver  and 
computer,  one  is  no  longer  available,  and  the  other  has  been  in  production  for  13  years,  a  long 
run  for  a  National  Instruments  product.  The  software  suite  that  runs  on  these  systems  is  a 
mixture  of  old  (FORTRAN)  and  new,  but  all  of  the  programs  operate  perfectly  well  on  the  latest 
computer  hardware. 

B2.  Description  and  Status  of  the  ITS  Receiver  Stations 
Gakona 

The  ITS  10  station  at  the  HAARP  site,  near  Gakona,  AK,  is  the  first  of  the  ITS  stations  to  be 
installed  as  part  of  the  HAARP  operation.  This  station  is  the  most  ideal  of  all  the  installations,  as 
RF  noise  is  minimal,  a  large  area  is  available  for  a  large  antenna  ground-plane,  and  the 
equipment  is  housed  in,  and  operated  from,  an  equipment  shelter  specifically  designed  for 
scientific  observation.  Besides  power,  the  facility  provides  telephone  service  and  high-speed 
Internet  connectivity.  The  on-call  availability  of  HAARP  staff  personnel  to  perform  occasional 
jobs  that  need  hands-on  operation  is  a  major  benefit  for  this  site  location. 

Delta 

The  Delta  station,  the  second  Alaska  ITS  station,  was  installed  in  December  1999  at  the 
Silver  Fox  Roadhouse.  Initially,  there  were  multiple  problems  with  receiver  tuning,  antenna 
cabling  and  very  noisy  telephone  service.  It  took  a  year  and  several  visits  to  resolve  the 
problems  but,  since  then,  the  Delta  station  has  been  very  reliable.  The  Silver  Fox  Roadhouse  is  a 
gas  station  /  convenience  store  with  rental  cottages,  mostly  used  by  recreational  hunters.  The 
receiver  is  located  in  a  small  back  room  away  from  activity  in  the  store.  The  antenna  is  mounted 
on  a  wooden  frame  with  a  mesh  ground-plane  similar  to  the  one  at  Gakona. 
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Delta  is  the  most  RF-quiet  location  of  all  the  sites.  Moreover,  since  a  fiber-optic  telephone 
cable  was  installed  along  the  highway  outside,  communication  has  been  very  reliable.  We  have 
the  option  of  installing  a  DSL  Modem  but,  at  present,  there  is  no  need  for  this.  For  the 
occasional  service  to  the  system,  we  have  always  had  the  assistance  of  the  owner  of  the  Silver 
Fox,  Mr.  Dan  Splain. 

Cordova 

The  Cordova  station  was  installed  in  November  of  2000.  It  was  first  installed  at  Cordova 
High  School  in  the  City  of  Cordova,  but  the  student-teacher  activities,  local  RF  interference,  and 
occultation  of  the  horizon  by  the  surrounding  mountains  led  to  the  abandonment  of  that  site.  The 
station  was  moved  to  Cordova  Auto  Rental,  a  small  business  operating  near  the  Cordova  airport, 
about  12  miles  from  the  town.  Except  for  a  nearby  aeronautical  navigation  beacon,  the  location 
is  very  radio-quiet.  (The  navigation  beacon  interference  is  removed  with  an  RF  filter.)  The 
horizon  is  low  in  all  directions  so  low-elevation  satellite  passes  are  obtainable.  The  antenna  is 
mounted  atop  an  8’  x  40’  steel  cargo  container,  which  provides  a  good  ground  plane.  The 
receiver  and  computer  are  inside  the  Rental  office,  inside  a  TV-enclosure  type  cabinet. 
Communication  with  the  Cordova  station  is  via  the  local  telephone  cooperative.  Cordova 
Telephone,  in  turn,  communicates  with  the  outside  world  entirely  by  satellite  communication. 
There  are  occasional  outages,  but  they  are  rarely  longer  than  a  few  hours. 

Other  receiver  outages  (not  caused  by  the  Telephone  system)  usually  require  the  help  of 
someone  on-site.  For  several  years,  the  owner  of  Cordova  Auto  provided  this  service.  However, 
after  he  sold  the  operation,  a  local  Cordova  resident,  Mr.  Jim  Yansant,  has  been  available  to  us  as 
the  need  arises. 

Arctic  Village 

The  station  at  Arctic  Village  was  installed  in  March,  2004.  It  began  as,  and  continues  to  be, 
a  difficult  environment  from  which  to  gather  data.  The  location  of  the  station  is  an  8’  x  8’ 
science  experiment  shelter  located  near  the  Arctic  Village  airport.  This  facility  is  owned  and 
operated  by  the  University  of  Alaska,  Fairbanks.  Indeed,  the  actual  installation  and  most 
subsequent  work  has  been  by  an  employee  of  UAF,  Ms.  Kay  Lawson. 

The  difficulty  at  Arctic  Village  is  twofold;  communication  with  the  outside  world,  and  radio 
interference.  The  facility  had  (and  has)  a  telephone  connection,  so  initially  it  was  believed  that 
communication  would  be  similar  to  Cordova.  That  is,  there  is  a  copper  telephone  connection 
from  the  receiver  station  to  the  village,  and,  then,  a  satellite  link  to  the  rest  of  the  world. 
However,  the  satellite  link  is  vastly  far  below  Cordova  in  terms  of  data  transmission  speed.  The 
solution  was  to  set  up  our  own  link  from  the  station  to  a  second  (Internet-only)  satellite  station 
within  the  village.  This  was  done  in  concert  with  an  Alaska-wide  program  to  provide  Internet 
access  to  all  Alaskans,  regardless  of  the  size  of  the  villages  they  lived  in.  This  system,  by  itself, 
works  adequately,  albeit  very  slowly  at  times.  The  problem  is  that  the  Internet  connection 
equipment  is,  itself,  a  source  of  interference  to  ITS  10  operations.  The  solution  here  was  to 
program  the  Internet  connection  equipment  to  power-down  during  a  satellite  pass.  This  mode  of 
operation  did  eliminate  the  equipment  interference,  but  sometime  during  the  setup,  testing,  and 
eventual  power-on  -  power-off  operating  mode,  another  source  of  interference  has  come  into 
being. 
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Arctic  Village  operation  remains  poor  in  terms  of  the  amount  of  useful  data  obtained  from 
the  site.  We  are  investigating  an  alternative  location  for  this  receiver,  perhaps  at  the  Toolik  Field 
Site,  operated  for  the  National  Science  Foundation  by  the  University  of  Alaska  at  Fairbanks. 
This  site  is  on  the  main  road  from  Fairbanks  to  Prudhoe  Bay,  and  would  provide  better 
communications  and  on-site  support. 

B3.  Receiver  Performance  /  Operational  Upgrades 

In  addition  to  serving  as  a  back-up,  there  has  always  been  an  effort  to  perform  various 
feature  upgrades  on  the  suite  of  receivers  without  specifically  bringing  them  in  for  service.  That 
is,  the  back-up  receiver  will  receive  whatever  upgrades  have  been  deemed  useful  and,  then, 
whichever  field  receiver  needs  service  will  be  replaced  with  the  upgraded  receiver.  So  far,  only 
once  has  an  otherwise  functional  receiver  been  replaced  solely  for  the  purpose  of  receiving  an 
upgrade.  The  computers,  while  old,  are  completely  adequate  to  their  task.  Software  revisions 
and  adjustments  have  been  made  and  installed  remotely  since  the  systems  were  installed. 

B4.  Routine  Operation  /  Equipment  Malfunction 

In  routine  operation,  the  only  interaction  required  is  the  weekly  sending  of  the  NORAD 
TLEs  to  all  stations.  Of  course,  with  four  stations,  all  operating  at  remote  locations,  the  out-of- 
routine  is  common.  The  usual  cause  of  interruption  of  the  data  stream  is  either  the  telephone 
connections  to  the  three  stations  that  depend  on  that  mode,  or  a  slowdown  in  the  data  archiving  at 
the  primary  server  computer  located  at  the  HAARP  facility.  This  kind  of  problem  will  often 
seem  to  fix  itself.  For  example,  the  local  telephone  company  will  fix  a  downed  line. 

Once  in  a  while,  an  individual  station  will  have  its  modem  become  inoperative,  and  a  power 
cycle  to  the  system  is  required.  This  is  the  usual  duty  of  the  on-site  people  that  we  depend  on. 
More  rarely,  there  will  be  a  component  failure  in  the  receiver  or  even,  more  rarely,  in  the 
computer.  When  this  happens,  the  disabled  equipment  is  shipped  to  NWRA's  laboratory  in 
Redmond.  If  the  difficulty  can  be  adequately  diagnosed  remotely,  the  replacement  receiver  / 
computer  will  first  be  shipped  to  the  problem  site,  and  the  on-site  person  will  remove  /  replace 
the  unit  and  return  the  defective  unit. 
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Appendix  C 
Classic  Riometer  Status 

The  following  report  on  operations  and  maintenance  on  the  HAARP  30  MHz  classic 
riometer  was  provided  by  NWRA  consultant  Mr.  Jens  Ostergaard. 

Cl.  Overview 

The  HAARP  30-MHz  classic  riometer  is  one  of  a  wide  range  of  diagnostic  instruments  at 
the  HAARP  research  site  in  Gakona  Alaska.  Such  riometers  have  been  a  standard  for 
measurement  of  ionospheric  absorption  for  many  years.  The  frequency  of  30  MHz  is  widely 
used,  so  data  from  different  stations  can  readily  be  compared.  The  HAARP  riometer  operates 
with  an  antenna  array  of  four,  5 -element  crossed  Yagi  antennas  pointed  at  Zenith.  The  3-dB 
beamwidth  is  approximately  60°  x  60°.  The  riometer  receiver  is  a  solid-state,  noise  balancing 
radiometer  receiver  with  a  digitally  controlled  noise  source  used  for  balancing. 

The  ionospheric  absorption  is  measured  as  the  difference  between  the  Galactic  noise 
background  and  the  measured  antenna  temperature.  The  Galactic  noise  background  and  its 
diurnal  variation,  the  Quiet  Day  Curve,  can  either  be  determined  from  registrations  during 
periods  where  the  ionosphere  is  very  quiet,  or  it  can  be  computed  from  maps  of  galactic  noise 
sources  and  computer  models  of  the  antennas  radiation  pattern.  Both  methods  have  been  used 
during  the  last  three  years.  Figure  Cl  shows  examples  of  calculated  Quiet  Day  Curves  for 
various  frequencies.  The  computed  Quiet  Day  Curve  is  very  accurate  during  the  mid-winter  and 
mid-summer  months,  where  the  physical  environment  around  the  antenna  is  dry  and  stable. 
However,  during  periods  of  frost  and  thaw,  especially  in  the  spring,  fairly  large  deviations  have 
been  registered  on  a  recurring  basis.  The  deviations  include  both  calibration  jumps,  and  changes 
in  the  form  of  the  Quiet  Day  Curve.  The  latter  suggests  changes  in  the  antenna  radiation  pattern 
with  changes  in  the  environment  around  the  antenna.  It  has  also  been  noted  that  hoar  frost  on  the 
antenna  elements  may  have  caused  some  of  the  deviations.  Figure  C2  illustrate  how  differences 
in  the  condition  of  the  soil  and  snow  accumulation  can  effect  the  calculated  Quiet  Day  Curve. 

Thus,  it  must  be  concluded  that  the  current  computational  model  of  the  Quiet  Day  Curve 
needs  refinement.  During  periods  of  instability,  Quiet  Day  Curves  derived  from  measurements 
on  ionospherically  quiet  days  have  been  used.  The  need  for  a  computed  Quiet  Day  Curve  arose 
during  the  peak  of  the  past  solar  cycle,  where  the  ionospheric  activity  was  so  high,  that  it  was 
not  possible  to  derive  a  Quiet  Day  curve  from  measurements  at  all.  The  past  three  years  have 
included  the  least  active  period  of  the  solar  cycle,  and  a  large  number  of  quiet  days  have  been 
available.  This  data  has  been  very  valuable  when  evaluating  the  absolute  accuracy  of  the 
riometer  receiver  and  the  antenna  system.  This  report  presents  some  of  the  results  of  the 
evaluation. 
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Figure  Cl.  Examples  of  Quiet  Day  Curves  calculated  for  the  HAARP  riometer  for  a  range  of  frequencies. 


Figure  C2.  Examples  of  the  impact  of  varying  soil  conditions  (left-side  plots)  and  varying  snow  accumulation 

(right-side  plots)  on  the  calculated  Quiet  Day  Curve. 

C2.  Operational  Routine 

The  riometer  operates  continuously,  and  the  Quiet  Day  Curve  is  included  in  the  operational 
software.  This  way,  the  derived  ionospheric  absorption  is  available  in  real  time.  Both  the 
measured  galactic  noise  data  and  the  Quiet  Day  Curve  are  stored,  and  the  absorption  can  be 
reworked  at  a  later  date,  if  the  Quiet  Day  Curve  is  deemed  too  inaccurate  for  a  particular  use. 
The  riometer  operation  is  checked  every  day,  and  calibration  jumps,  as  well  as  adjustments  to  the 
Quiet  Day  Curve  are  performed  remotely  as  needed.  The  best  calibration  check  is  a  comparison 
of  the  Quiet  Day  Curve  with  measurements  during  quiet  days.  Quiet  Day  Curve  changes  are  not 
introduced  until  the  measurements  have  been  observed  for  some  days,  and  it  has  become  evident 
that  the  deviations  are  caused  by  the  instrument,  not  ionospheric  activity.  Weak  absorption 
events,  not  associated  with  magnetic  disturbances,  can  be  observed  at  HAARP  during  late 
afternoon  hours  and  early  evening  hours.  These  events  are  caused  by  particles  dumped  out  of  the 
radiation  belts.  (Dr.  Lee  Snyder,  private  communication).  It  is  important  not  to  interpret  such 
events  as  changes  in  the  Quiet  Day  Curve. 

A  site  visit  is  performed  once  per  year.  The  maintenance  work  includes  antenna  leveling. 
This  is  necessary,  as  the  guy  anchors  embedded  in  the  permafrost  move  slightly  every  year.  The 


64 


return  loss,  i.e.  SI  1  are  measured  in  the  frequency  band  29.7  MHz  -  31  MHz  for  all  antennas  and 
feed  systems.  A  high  return  loss  is  the  most  suitable  way  to  determine  if  bad  connections  or 
corrosion  exist  in  the  antennas  or  the  feed  system.  Deficiencies  are  repaired  as  necessary. 
Finally,  the  riometer  receiver  is  calibrated  with  a  variable  noise  source. 

C3.  Mechanical  Failures  and  Problems 

The  climate  at  HAARP  is  very  hard  on  equipment  located  outdoors.  The  riometer  antenna 
array  has  experienced  a  number  of  mechanical  problems  over  the  years.  The  antenna  elements, 
as  such,  have  held  up  well,  and  they  should  last  for  many  more  years.  There  have,  however,  been 
a  number  of  failures  with  the  baluns  in  the  feed  system.  The  antennas  require  a  1 :4  step-up  balun 
between  the  feeder  cable  and  the  fed  element.  The  original  baluns  were  made  with  very  thin 
coaxial  cable  potted  in  a  hard  epoxy.  These  baluns  did  have  a  number  of  failures  most  likely  due 
to  the  internal  coax  cable  being  so  constrained  that  the  thin  inner  conductor  was  tom  by  the 
extreme  temperature  variations  at  the  site.  The  baluns  were  replaced  a  few  years  ago  with 
standard  half-wave  coaxial  baluns  made  from  RG213  coaxial  cable  connected  with  N-connectors 
to  a  watertight  junction  box.  All  junctions  are  tightened  hard  and  covered  with  shrink  tubing  The 
new  baluns  have  lasted  well  for  the  past  three  years.  The  baluns  need  maintenance  during  the  site 
visits.  The  N  connectors  can  work  loose,  and  the  shrink  tubing  does  deteriorate  with  time  and 
needs  to  be  replaced.  The  coax  cables  leading  from  the  antennas  to  the  riometer  receiver  lasted 
10  years,  but  had  to  be  replaced  in  the  spring  of  2008  due  to  cracks  in  the  cable  jacket. 

It  was  found  in  the  spring  of  2008  that  Pad  3  does  not  drain  well  during  the  spring  thaw. 
Large  puddles  of  water  covered  the  pad  for  a  period  of  many  weeks.  The  pad  surface  has  become 
uneven  and  water  tends  to  congregate  in  the  low  spots.  One  such  spot  is  located  under  the  box 
protecting  the  riometer  receiver,  and  1/2”  of  water  was  found  inside  the  box  at  the  end  of  the 
spring.  The  low  spot  will  need  to  be  filled  in,  and  the  protective  box  placed  on  a  bed  of  cinder 
blocks  to  keep  it  above  the  water  during  the  thaw.  All  cable  transitions  into  the  box  have  been 
checked  and  resealed. 

The  moisture  in  the  box,  in  combination  with  the  puddles  on  the  pad,  caused  problems  with 
the  power  cables  leading  from  the  equipment  shelter  to  the  riometer  antennas.  The  GFI  relay 
protecting  the  cable  dropped-out  daily.  The  cable  was  inspected  and  tested  for  galvanic  leakage. 
The  cable  is  in  good  condition,  and  no  leakage  was  found.  The  GFI  relay  has  been  replaced.  The 
GFI  relay  dropouts  stopped  when  the  pad  dried  out  in  the  late  spring. 

C4.  Evaluation  of  Accuracy 

The  resolution  ( R )  of  a  noise  balancing  radiometer  receiver  is  calculated  as: 


where  T  is  the  total  noise  figure,  B  is  the  bandwidth,  and  tt  is  the  integration  time.  The  total  noise 
figure  for  the  HAARP  riometer  including  feeder  losses  is  5  dB  or  approximately  600°K,  and  the 
bandwidth  is  30  kHz.  At  the  maximum  recorded  noise  level  of  22,000°K,  the  resolution  is  better 
than  0.1  dB.  Similarly,  at  13-dB  absorption,  the  resolution  is  0.1  dB;  both  for  an  integration  time 
of  1  second.  A  better  resolution  is  obtained  with  longer  integration  times.  The  long-term, 
absolute  accuracy  with  the  current  riometer  receiver  is  not  as  good.  The  influence  of  calibration 
jumps  and  presumed  antenna  pattern  changes  have,  during  some  periods,  introduced  a  0.2-dB  - 
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0.4-dB  difference  between  the  Quiet  Day  Curve  and  measurements  during  quiet  ionospheric 
conditions.  This  is  not  satisfactory.  The  goal  for  the  overall  accuracy  of  the  HAARP  riometer  is 
0.1  dB.  The  following  paragraphs  examine  some  of  the  sources  of  these  variations. 

C5.  Properties  of  the  Antenna  Array 

Extensive  modeling  of  the  antenna  array  radiation  pattern  has  been  performed  with  NEC. 
The  antennas  useful  frequency  range  has  been  determined  to  cover  29.7  MHz  -  31  MHz.  The 
gain  in  this  range  varies  0.3  dB,  and  the  low-angle  side  lobes  are  well  suppressed.  The  main  lobe 
deteriorates  rapidly  with  an  increase  in  frequency  above  3 1  MHz.  At  32  MHz,  the  front  lobe  and 
the  low-angle  side  lobe  have  equal  magnitude  and  the  Quiet  Day  Curve  has  degenerated.  The 
gain  figures/beamwidth  calculations  have  been  corroborated  by  operating  the  riometer  in  a 
frequency  sweeping  mode. 

A  free  space,  NEC  model  at  30  MHz  shows  that  the  radiation  pattern  has  a  respectable  ffont- 
to-back  ratio  of  20  dB,  and  a  front  gain  of  14  dBi.  The  main  beam  width  is  slightly  larger  than  60 
deg.  Models,  including  a  ground  under  the  antenna,  show  a  gain  variation  from  14.4  dBi  to  15 
dBi  dependent  on  the  soil  properties  used.  The  highest  gain  is  associated  with  a  perfect  ground, 
the  lowest  with  a  dry,  sandy  soil  model.  Corresponding  small  variations  in  the  beam  width  can 
also  be  seen.  Quiet  Day  Curves  have  been  computed  based  on  these  models,  and  a  small 
deviation  can  be  seen  at  the  Galactic  maximum  part  of  the  curve.  The  low  antenna  temperature 
range  of  the  Quiet  Day  Curves,  where  the  antenna  points  farthest  away  from  the  ecliptic  plane, 
show  a  larger  variation.  The  noise  distribution  over  the  antenna  aperture  is  very  uniform  in  this 
area.  A  slight  variation  in  the  gain  does  not  change  the  antenna  temperature,  but  the  associated 
widening  of  the  beam  does.  A  larger  effect  is  seen  if  the  distance  between  the  ground  plane  and 
the  antenna  is  decreased.  Again,  a  perfect  ground,  or  a  wet  snow  pack  during  the  thaw  period 
has  a  larger  effect  than  a  dry  and  frozen  snow/soil  model.  This  could  explain  some  of  the 
noticeable  variations  in  the  form  of  the  Quiet  Day  Curve  during  the  thaw  period.  Further  analysis 
will  be  needed  to  confirm  this  hypothesis.  It  is  tempting  to  install  a  ground  screen  under  the 
antenna.  The  materials  are  available  at  the  site.  However,  such  a  screen  would  be  located  under 
the  winter  snow  accumulation,  and  would  not  remove  the  influence  of  a  wet  snow  pack.  Further 
modeling  will  be  necessary  to  examine  the  effects  of  a  combination  of  a  ground  screen  and  the 
snow  cover. 

C6.  Properties  of  the  Riometer  Receiver 

The  advantage  of  the  Noise  Balancing  Principle  radiometer  receiver  is  that  gain  variation 
and  variation  in  the  receiver’s  noise  figure  are  common  to  both  the  antenna  branch  and  the  noise 
generator  branch.  This  eliminates  calibration  concerns  for  the  receiver  circuits,  detectors  and 
filters.  It  is  still  necessary  to  account  for  the  loss  in  the  antenna  feeder,  but  that  is  easily 
measured,  and  supposedly  very  stable,  unless  the  cable  is  damaged.  Other  sources  of  calibration 
drift  and  instability  exist,  however.  The  solid-state  noise  source  used  for  noise  balancing  is  an 
avalanche-type  diode  with  an  Excess  Noise  Ratio  drift  of  0.01  dB/C.  That  is  very  little,  but  over 
an  ambient  temperature  range  of  +30°C  to  -40°C,  this  could  amount  to  almost  0.1  dB.  The 
HAARP  riometer  is  enclosed  in  a  thermostatic  box  that  eliminates  this  drift.  Other  sources  are 
offset  voltages  in  the  feed  back  loop  integrator  that  vary  with  temperature  and  time,  and  gain 
instabilities  in  the  receivers  IF  amplifier.  Riometers  are  commonly  designed  as  single  conversion 
receivers  with  a  10.7-MHz  IF  frequency.  To  obtain  front-end  robustness,  most  of  the  signal 
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amplification  is  placed  in  the  IF  amplifier.  This  is  also  the  case  for  the  HAARP  riometer,  where 
the  robustness  is  needed  due  to  the  proximity  of  the  transmitter  array.  A  dual  conversion 
receiver,  where  more  of  the  signal  amplification  is  moved  to  either  baseband  or  a  low  frequency 
RF  amplifier,  could  eliminate  risk  of  IF  amplifier  instabilities.  A  lower  IF  amplifier  gain  and  a 
more  sophisticated  baseband  amplifier/detector/integrator  is  part  of  the  riometer  receiver  upgrade 
currently  underway. 

The  riometer  data  user  community  has  expressed  interest  in  a  better  time  resolution  than  the 
current  one  minute.  A  sample  of  data  with  a  time  resolution  of  one  second  has  been  collected. 
Examples  of  both  quiet  and  highly  disturbed  periods  are  included  in  the  data.  Analysis  of  this 
data  is  still  pending. 

Cl.  Conclusion 

Operational  experiences  with  a  new  riometer  receiver  have  been  accumulated  throughout  the 
last  three  years.  This  period  includes  the  low  sunspot  part  of  the  solar  cycle,  and  extended 
periods  of  quiet  ionospheric  conditions  have  been  very  useful  in  evaluating  and  improving  the 
overall  accuracy  of  the  HAARP  classic  30-MHz  riometer.  It  has  been  proven  that  the  use  of  a 
computed  Quiet  Day  Curve  is  feasible.  However,  seasonal  calibration  jumps  still  reduce  the 
absolute  accuracy  more  than  desirable.  The  report  has  outlined  the  daily  and  annual  operation 
and  maintenance  routines;  and  has  presented  a  number  of  factors  influencing  the  measurement 
accuracy.  The  improvement  work  is  ongoing,  and  will  include  work  on  the  ground  screen,  the 
antenna  supports  and  revisions  of  the  riometer  receiver. 
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Appendix  D 

Ionospheric  Scintillation  Modeling  Work 


This  task  was  devoted  to  two  major  efforts.  The  first  effort  involved  the  use  of 
measurements  of  total  electron  content  (TEC)  from  ALTAIR  to  construct  a  model  of  the 
ionosphere  for  application  to  problems  involving  space  based  radar.  The  second  effort  involved 
the  application  and  validation  of  a  Matlab  version  of  the  multiple  phase  screen  (MPS) 
propagation  code. 

Dl.  Model  of  the  Ionosphere  Based  on  Radar  Data 

Consider  the  1988  PEAK  experiment.  The  1988  PEAK  (Propagation  Effects  Assessment  - 
Kwajalein)  experiment  collected  ALTAIR  dual  frequency  (VHF  and  UHF)  pulse-by-pulse  range 
measurements  during  tracks  of  calibration  spheres.  The  total  electron  content  during  these  tracks 
was  obtained  from  the  range  measurements  using  the  difference  in  the  ranges  measured  at  VHF 
and  UHF,  [Knepp  and  Houpis,  1992],  In  1988,  the  radar  I&Q  data  was  also  analyzed  to  obtain 
values  of  the  scintillation  index  and  decorrelation  time  during  track  of  COSMOS  1427,  a  large 
Russian  calibration  sphere.  For  this  track  (referred  to  as  PEAK2G),  we  recently  found  a  binary 
file  (written  in  1989)  with  machine  precision  TEC  measurements.  The  accuracy  of  these 
measurements  is  limited  by  the  radar  receiver  noise,  not  by  the  computer  word  length. 

In  the  following  we  describe  a  straightforward  technique  to  utilize  this  dual-frequency 
ALTAIR  radar  track  data  to  construct  an  ionospheric  model.  Given  the  model,  we  then  apply  an 
MPS  code  to  predict  the  scintillation  index  and  decorrelation  time.  These  model  predictions  are 
then  compared  to  the  actual  measurements  of  scintillation  index  and  decorrelation  time.  This  and 
other  similar  comparisons  are  useful  to  validate  the  model. 

The  PEAK  experiment  that  was  conducted  in  1988  was  originally  planned  to  be  the  first  of 
two  data  collections,  with  the  second  to  occur  two  years  later  at  solar  maximum.  The  Cold  War 
ended  in  1989,  and  the  second  data  collection  experiment  was  canceled  by  the  Defense  Nuclear 
Agency.  Mission  Research  Corporation  was  directed  to  get  as  much  as  possible  out  of  the  1988 
data  and  quickly  wrap  up  the  program.  For  this  reason,  the  satellite  track  referred  to  as  PEAK2G 
has  the  only  set  of  measurements  of  scintillation  index,  decorrelation  time,  and  coherence 
bandwidth.  These  measurements  are  used  in  the  comparison  discussed  here. 

In  the  processing  described  below,  the  ALTAIR  TEC  data  is  used  to  construct  a  model  of 
the  ionosphere.  Then,  the  MPS  code  uses  this  model  to  calculate  the  propagation  conditions 
observed  at  ALTAIR  using  the  measured  TEC  as  the  basis  of  the  MPS  phase  screen.  Finally,  the 
experimentally  measured  scintillation  index  and  the  decorrelation  time  are  compared  to  those 
values  calculated  from  the  MPS  simulation.  This  comparison  allows  one  to  determine  the 
usefulness  and  accuracy  of  the  process  used  to  generate  the  phase  screen  representation  of  the 
ionosphere. 

Fortunately,  a  file  with  some  relevant  TEC  data  from  PEAK2G  was  available  from  a  backup 
CD.  This  file  is  analyzed  here.  The  file  is  a  Vax  binary  file  and  contains  TEC  measurements 
and  sufficient  radar  information  to  generate  radar  range,  azimuth,  and  elevation.  Matlab  is  used 
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to  quickly  develop  code  to  process  this  data.  The  data  is  described  in  several  reports  and  papers 
including  Knepp  and  Houpis  [1991  and  1992]  and  Houpis  and  Nickisch  [1991], 

The  data  analyzed  here  consist  of  the  TEC  measured  pulse  by  pulse  from  simultaneous  VHF 
and  UHF  range  measurements  obtained  during  the  track  of  a  Russian  calibration  sphere  in  a 
polar  orbit  at  an  altitude  of  420  km.  The  goal  of  this  work  is  to  use  the  TEC  measured  in  this 
fashion  to  generate  a  realistic  phase  screen  representation  of  the  ionosphere.  This  phase  screen  is 
then  used  in  an  MPS  propagation  simulation  to  “recreate”  the  scintillation  measurements 
observed  by  ALTAIR.  In  the  following,  we  compare  the  actual  data  (measured  scintillation 
index  and  decorrelation  time)  with  that  produced  by  this  phase  screen  model.  Once  validated, 
such  a  process  could  then  be  applied  to  other  transmission  frequencies  and  geometries  such  as 
used  by  space-based  synthetic  aperture  radar. 

The  presence  of  receiver  noise  must  be  taken  into  account.  Although  the  signal-to-noise 
ratio  varies  during  the  track,  typical  values  are  20  dB  at  VHF  and  30  dB  at  UHF.  At  20-dB  VHF 
SNR  and  30-dB  UHF  SNR,  the  range  standard  deviations  (due  to  noise)  at  VHF  (7-MHz 
bandwidth)  and  UHF  (18-MHz  bandwidth)  are  1.52  and  0.18  m,  respectively.  Using  the  range 
measurements  to  obtain  the  TEC,  the  TEC  standard  deviation  is  0. 1 1  TEC  units. 


Converting  the  TEC  directly  into  phase  for  a  screen  at  UHF  (422  MHz)  gives  a  standard 
deviation  for  the  phase  of  2.2  radians  or  126  degrees.  At  VHF,  the  standard  deviation  is  higher 
by  a  factor  of  (422/155).  This  large  an  error  cannot  be  handled  correctly  in  the  MPS  propagation 
simulation  which  requires  smoothly  varying  phase  from  sample  to  sample.  (The  criteria  for 
application  of  numerical  MPS  techniques  is  discussed  in  Knepp  [1983].)  Thus  some  creative 
processing  of  the  original  TEC  data  must  be  developed  to  obtain  smooth  values  of  phase. 


File:  2GVTEC.SEQ  Case:  19  Pass:  Geometry 
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Figure  Dl.  Geometry  of  pass  PEAK2G. 


Figure  Dl  shows  the  trajectory  for  PEAK2G.  The  red  x  shows  the  location  of  the  ALTAIR 
radar,  the  blue  curve  shows  the  latitude  and  longitude  of  the  satellite,  and  the  green  curve  shows 
the  latitude  and  longitude  of  the  ionospheric  intercept  point  at  the  assumed  altitude  of  350  km. 
In  order  to  simulate  this  geometry  with  a  phase  screen,  the  data  here  is  transformed  to  obtain  the 
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TEC  as  a  function  of  distance  along  the  trajectory  and  perpendicular  to  the  direction  of  the 
earth’s  magnetic  field.  This  is  accomplished  by  converting  the  latitude  and  longitude  of  the 
trajectory  to  magnetic  latitude  and  longitude  using  an  earth-centered  dipole  magnetic  field 
model.  Then,  the  local  angle  with  respect  to  magnetic  North  is  obtained  along  the  trajectory. 
For  the  trajectory  in  PEAK2G,  this  angle  is  approximately  constant  during  the  entire  track,  so  the 
TEC  data  is  easily  obtained  as  a  function  of  distance  perpendicular  to  the  earth’s  magnetic  field. 
The  data  consisting  of  TEC  versus  distance  is  the  basis  of  the  single  phase  screen  representation 
of  the  ionosphere  to  be  used  here. 

For  PEAK2G,  the  angle  of  the  trajectory  with  respect  to  magnetic  north  is  26.5  degrees.  To 
obtain  the  distance  along  the  trajectory  perpendicular  to  magnetic  north,  one  simply  divides  the 
measured  distance  along  the  trajectory  by  the  sine  of  26.5  degrees.  The  temporal  spacing 
between  radar  data  samples  is  variable.  Consequently,  the  spatial  separation  is  also  variable. 
The  average  spatial  separation  for  PEAK2G  is  41.2  meters  and  the  maximum  separation  is  261.3 
meters.  Since  the  separations  are  not  constant,  interpolation  must  be  applied  in  order  to  obtain 
the  constant  spacing  necessary  to  the  FFT  algorithm  in  the  MPS  code.  Of  course,  some  filtering 
of  the  noise  is  required  prior  to  the  interpolation.  Otherwise,  the  interpolation  process  will  tend 
to  increase  the  contribution  of  the  noise. 

The  other  required  quantity  for  the  MPS  code  is  the  propagation  distance  for  the  effective 
phase  screen.  The  MPS  code  presently  models  only  plane  wave  propagation  a  distance  Z  beyond 
a  phase  screen.  By  the  simple  lens  law,  spherical  propagation  is  well-modeled  using  the 
expression  Z*  =  Z1Z2X  Zi+  Z2)  where  Zi  is  the  distance  from  the  transmitter  to  the  phase  screen 
and  Z2  is  the  distance  from  the  phase  screen  to  the  receiver. 

The  MPS  code  was  initially  intended  for  one-way  propagation  from  a  satellite  to  a  ground 
receiver.  To  model  two-way  radar  propagation,  the  calculated  voltage  is  squared  (Eq.  12  in 
Knepp  and  Houpis,  [1991]). 

For  the  PEAK2G  data,  Figure  D2  shows  the  values  of  the  distances  as  a  function  of 
ALT  AIR  radar  time  during  the  track.  The  value  of  Z*  ranges  from  137  to  187  km  during  the 
track.  In  the  MPS  calculation,  the  mean  value  of  158  km,  obtained  as  the  average  of  all  the 
individual  values  during  the  trajectory,  is  used  to  specify  the  constant  distance  from  the  phase 
screen  to  the  receiver  plane. 

The  measured  value  of  slant  TEC  is  multiplied  by  the  classical  electron  radius  and  the 
wavelength  to  calculate  the  phase  at  the  radar  transmission  frequency.  Then,  to  remove  some  of 
the  noise,  we  smooth  this  phase  by  replacing  the  values  with  the  output  of  a  5-point  running 
average.  On  average,  this  fdters  the  data  over  a  distance  of  5  times  41.2  meters,  or  206  meters, 
which  is  comparable  to  the  Fresnel  zone.  The  Fresnel  zone  is  the  square  root  of  the  wavelength 
times  Z*,  which  at  VHF  (155  MHz)  is  552  meters.  Irregularities  with  sizes  on  the  order  of  the 
Fresnel  zone  contribute  most  to  scintillation.  Such  filtering  will  tend  to  remove  structure  that  is 
on  the  order  of,  and  smaller  than,  the  filter  length.  Too  much  filtering  of  the  phase  will  thus  tend 
to  remove  scintillation-producing  irregularities.  To  better  understand  the  effect  of  this  filtering, 
we  performed  a  numerical  experiment  using  the  MPS  simulation. 
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File:  2GVTEC.SEQ  Case:  19  PhaseScreen:  Distance 


Figure  D2.  Propagation  distances  during  PEAK2G. 


The  effect  of  this  filtering  depends  on  the  spectrum  of  the  phase.  For  the  phase  spectrum  of 
PEAK2G  which  has  a  slope  of  roughly  2.5,  we  found,  by  numerical  experimentation  (discussed 
below),  that  the  effect  of  a  5-point  smoothing  filter  is  small  for  the  VHF  case.  For  the  UHF  case, 
with  its  smaller  Fresnel  zone  of  335  meters,  the  effect  of  the  5-point  filter  is  more  important  but 
still  minor. 

To  continue  the  description  of  the  processing,  the  filtered  phase  versus  distance  is  then 
linearly  interpolated  to  obtain  a  series  of  equally  spaced  phase  values  separated  by  a  distance  of 
20.6  meters.  For  PEAK2G,  the  original  TEC  data  comprised  a  series  of  31,087  values.  The 
interpolation  step  described  above  increases  the  number  of  values  by  about  a  factor  of  2  to 
62,175  values. 

The  final  step  in  the  creation  of  a  phase  screen  model  is  to  increase  the  total  number  of 
points  to  a  power  of  2  (131,072  points).  For  the  phase  in  the  additional  points,  we  simply  used  a 
linear  (straight-line)  value  connecting  the  last  point  to  the  first  point  in  the  periodic  repetition  of 
the  phase  screen.  This  linear  function  connects  the  first  and  last  original  phase  values  so  that 
there  is  no  discontinuity  in  periodic  repeats  of  the  phase  screen.  Figure  D3  shows  the  complete 
phase  screen  at  VHF  that  is  used  to  generate  the  propagation  results  presented  in  this  report. 

Figure  D4  gives  a  comparison  of  the  VHF  phase  power  spectrum  after  the  filtering  and 
interpolation  steps  described  above.  The  right-hand  frame  shows  the  case  of  no  filtering.  The 
left-hand  frame  shows  the  effect  on  the  phase  PSD  of  the  5-point  running  average  filter.  The 
appended  linear  portion  of  the  phase  screen  is  not  used  to  calculate  this  power  spectrum.  The 
major  differences  appear  at  values  of  wavenumber,  q,  of  0.004  m_1and  above,  and  are  caused  by 
the  5 -point  running  average  filter. 
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Figure  D4.  Effect  of  filtering  on  the  phase  PSD. 


Figure  D5  shows  the  results  of  a  number  of  statistical  MPS  calculations  of  the  scintillation 
index  and  the  decorrelation  distance  for  a  simulation  of  the  PEAK2G  data.  These  results  are 
generated  with  a  number  of  statistically  independent  phase  screens  designed  to  simulate  the 
power  spectrum  observed  during  PEAK2G.  The  phase  screen  shown  in  Figure  D3  is  not  used  to 
generate  the  data  in  Figure  D5.  For  this  simulation,  a  phase  power  spectrum  with  a  slope  of  -2.5 
is  used  to  match  the  measured  power  spectrum.  The  propagation  distance,  frequency,  grid  size 
and  grid  spacing,  and  other  MPS  input  parameters  match  the  PEAK2G  measurements. 
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Effect  of  running-average  filter 


Figure  D5.  MPS  simulation  of  the  effect  of  the  filter  on  the  measured  scintillation  index  and 
decorrelation  time. 

Figure  D5  shows  the  effects  of  filters  of  various  lengths  that  are  used  to  modify  the 
statistical  MPS  phase  screen  prior  to  the  propagation  calculation.  The  quantities  plotted  are  the 
scintillation  index  divided  by  the  scintillation  index  for  the  case  of  no  phase  filtering  and  the 
decorrelation  distance  divided  by  the  decorrelation  distance  for  the  case  of  no  filtering.  These 
two  quantities  are  plotted  versus  the  filter  length.  The  plotted  parameters  increase  with 
increasing  filter  length  because  the  filter  is  removing  the  smaller  scale  irregularities  that  produce 
scintillation  but  leaving  larger  scale  structure  that  produces  large  lenses  that  give  strong  focusing 
and,  hence,  very  large  enhancements  in  power.  These  enhancements  increase  the  value  of 
scintillation  index.  As  the  filter  length  increases,  the  overall  phase  variance  decreases,  which 
reduces  the  total  scintillation  and  thus  increases  the  value  of  the  decorrelation  distance. 

The  results  presented  in  Figure  D5  show  that  the  filtering  does  indeed  impact  the  values  of 
the  scintillation  index  and  decorrelation  distance  and  should  be  avoided,  if  possible. 

In  the  following,  we  present  a  comparison  of  the  ALTAIR  data  (for  RCS,  scintillation  index, 
and  decorrelation  time)  versus  the  output  of  the  MPS  propagation  calculation  that  uses  the 
filtered  TEC  data  shown  in  Figure  D3  as  the  phase  screen  representation  of  the  ionosphere.  A 
comparison  of  Figures  D6  and  D7  show  the  power  versus  time  of  the  measured  ALTAIR  data 
and  the  MPS  simulation.  For  these  figures  alone,  a  small  amount  (20-dB  SNR)  of  Gaussian 
noise  is  added  to  the  output  electric  field  of  the  MPS  calculation.  This  additive  noise,  for 
plotting  purposes  only,  reduces  the  depth  of  some  of  the  deepest  fades  and  gives  a  noisy  power  in 
the  portion  of  the  MPS  grid  where  the  phase  screen  is  linear  (1200  to  3000  km  in  Figure  D3). 
The  level  of  scintillation  is  similar  in  the  MPS  simulation  and  the  actual  data.  The  major 
difference  in  the  figures  is  that  the  period  in  the  actual  data  from  about  550  to  650  seconds  is 
relatively  scintillation-free,  while  the  MPS  simulation  shows  severe  scintillation  during  that  same 
period.  This  indicates  that  the  5-point  moving  average  filter  did  not  remove  all  the  phase  (or 
TEC)  noise.  As  discussed  above,  additional  filtering  (more  than  5  points)  would  also  remove  the 
important  phase  structure  at  the  Fresnel  zone  size,  which  gives  the  most  important  contribution 
to  amplitude  scintillation. 
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Figure  D6.  RCS  measured  at  VHF  during  PEAK2G. 
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Figure  D7.  MPS  calculation  of  VHF  power  and  phase. 

Figures  D8  and  D9  show  the  value  of  scintillation  index  versus  radar  time  from  the  actual 
data  and  from  the  MPS  simulation.  Roughly  the  same  number  of  points  is  used  for  both 
calculations.  The  overall  appearance  of  the  curves  is  similar  with  the  same  degree  of  scintillation 
severity.  However,  the  simulation  fails  to  predict  the  reduction  in  scintillation  starting  around 
550  seconds.  In  the  original  data  processing  used  to  produce  Figure  D8,  1000  points  are  used  to 
calculate  a  single  value  of  scintillation  index.  This  1000-point  window  is  then  moved  in  time, 
and  additional  values  of  scintillation  index  are  calculated. 
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Figure  D8.  Scintillation  index  measured  at  VHF  during  PEAK2G. 

S4  Case:  55  L:  PEAK2G  Filt:  5  MHz:  155.5 


Figure  D9.  MPS  calculation  of  scintillation  index  at  VHF. 

In  the  simulation,  2000  points  are  used  for  this  same  process  because  the  effective  temporal 
spacing  in  the  simulation  is  twice  that  of  the  original  ALTAIR  data.  The  value  of  scintillation 
index  obtained  from  these  points  is  plotted  in  the  figure  at  the  center  of  this  data  window.  The 
window  is  moved  slightly  (one-quarter  of  its  length)  and  the  process  is  repeated.  Thus  the  values 
of  scintillation  index  shown  in  the  figure  are  correlated  and  the  two  figures  are  directly 
comparable. 

Figures  DIO  and  Dll  show  the  value  of  the  decorrelation  distance  measured  using  1024- 
point  FFTs  for  the  data,  and  2048-point  FFTs  for  the  simulation.  A  comparison  of  these  two 
figures  indicates  that  the  maximum  level  of  scintillation  is  similar  for  the  data  and  the  simulation, 
characterized  by  a  minimum  value  of  decorrelation  time  of  0.01  sec  in  both  cases.  The  curves 
look  similar,  but  the  value  of  decorrelation  time  in  the  data  is  increasing  later  in  the  pass  at  a 
much  higher  rate  than  in  the  simulation. 
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Figure  DIO.  Decorrelation  time  measured  at  VHF  during  PEAK2G. 
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Figure  Dll.  MPS  calculation  of  decorrelation  time  at  VHF. 


The  phase  screen  is  scaled  to  the  UHF  frequency  of  422  MHz  and  the  propagation 
calculation  proceeds  exactly  as  described  above.  Figure  D12  shows  the  power  at  UHF  from  the 
actual  ALT  AIR  track  and  Figure  D13  shows  the  power  measured  in  the  MPS  simulation.  A 
comparison  of  the  two  figures  shows  a  basic  similarity  of  results. 
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Figure  D12.  RCS  at  UHF  measured  during  PEAK2G. 
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The  overall  results  indicate  that  this  type  of  analysis  could  be  a  useful  approach  in  applying 
the  ALTAIR  data  to  other  radar  problems.  However,  the  contribution  to  noise  in  the  TEC 
measurements  remains  important  and  cannot  be  removed  or  filtered  very  well. 

The  effects  of  the  noise  could  be  overcome  by  another  approach.  ALTAIR  also  collects 
coherent  amplitude  and  phase  data.  It  is  this  data  that  was  used  to  measure  the  decorrelation  time 
shown  above.  Consider  the  phase  history  that  consists  of  the  measured  phase  as  a  function  of 
target  position  during  a  satellite  pass  observed  by  the  radar.  There  are  two  contributions  to  the 
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phase;  one  due  to  the  target  range  and  the  second  due  to  the  ionization.  The  contribution  due  to 
range  could  be  accurately  estimated  by  using  a  trajectory  fitting  code  to  filter  the  range  data  and 
obtain  the  phase  history  due  solely  to  range  effects.  Subtraction  of  the  phase  due  to  range  from 
the  measured  phase  would  yield  the  phase  due  to  ionization.  This  remainder  could  be  used  to 
generate  accurate  phase  screen  realizations  of  the  ionosphere. 

Similar  techniques  have  been  tried  using  differential  GPS  phase  and  in  the  analysis  of 
Transit  satellite  data  to  measure  TEC.  The  phase  obtained  in  this  fashion  is  a  good  measure  of 
the  ionospheric  phase  only  if  the  scintillation  index  of  the  data  is  less  than  0.75  and  the  sampling 
rate  is  sufficiently  high  [Knepp,  2004],  The  sampling  rate  criteria  appear  to  be  met  by  the 
ALT  AIR  data. 

D.2  Example  and  Validation  of  the  MPS  Code 

In  this  section  we  describe  the  application  of  the  Multiple  Phase  Screen  (MPS)  code  to  wide 
bandwidth  calculations  of  propagating  signals.  This  report  presents  results  from  a  two- 
dimensional  MPS  calculation  that  computes  numerical  realizations  or  sample  functions  of  the 
electric  field  after  propagation  through  the  ionosphere.  The  MPS  code  solves  the  parabolic  wave 
equation  and  allows  for  direct  computation  of  realizations  of  the  received  signal.  The  ionization 
is  represented  by  a  series  of  random  phase  screens  that  characterize  the  severity  and  spectrum  of 
the  electron  density  fluctuations.  The  user  can  also  design  their  specific  phase  screen  to 
represent  a  deterministic  ionosphere.  For  wide  bandwidth  signals,  the  MPS  code  is  exercised  for 
many  frequencies  over  the  bandwidth  of  the  propagating  signal;  Fourier  transform  techniques  are 
used  to  obtain  the  propagating  signal  in  the  time  domain. 

The  MPS  simulation  is  quite  general,  and  may  be  applied  to  problems  involving  numerous, 
separated,  layers  of  ionization  characterized  by  spatially  varying  electron  density  power  spectra. 
MPS  techniques  can  handle  all  levels  of  ionospheric  disturbances  from  the  least  severe,  where 
only  small  phase  fluctuations  occur,  to  the  most  severe  case  of  frequency  selective  scintillation. 

To  validate  the  MPS  simulation  in  this  report,  a  number  of  calculations  are  repeated  and  the 
realizations  are  averaged  to  obtain  the  two-frequency,  two-position  mutual  coherence  function, 
which  is  the  Fourier  transform  of  the  generalized  power  spectrum.  The  MPS  results  are 
compared  with  analytic  results  obtained  using  Nickisch’s  [1992]  phase  screen  diffraction  method 
(PDM). 

Small  Scale  Ionization 

Small-scale  electron  density  irregularities  can  cause  random  variations  of  the  amplitude, 
phase,  angle-of-arrival,  and  time-of-arrival  of  a  propagating  signal.  These  variations  are 
collectively  referred  to  as  scintillation.  Their  impact  on  a  radar  or  communications  system  is  a 
function  of  the  in-situ  ionospheric  structure,  the  characteristics  of  the  propagating  signal,  and  the 
radar  or  communications  system  design.  Since  the  index-of-refraction  in  a  plasma  is  frequency- 
dependent,  signals  at  lower  transmission  frequencies  and  wider  bandwidth  will  generally 
experience  more  severe  effects  in  comparison  to  signals  with  higher  transmission  frequencies. 

The  severity  of  the  ionospheric  structure  depends  on  the  local  ionospheric  weather 
conditions.  In  the  equatorial  region,  strong  scintillation  can  be  important  at  VHF  and  UHF. 
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Equatorial  scintillation  begins  in  post-sunset  and  can  last  past  local  midnight.  In  the  polar 
regions,  scintillation  can  be  important  anytime  of  the  day  or  night. 

MPS  Calculation 

In  this  report,  we  apply  a  multiple  phase  screen  (MPS)  technique  [Knepp,  1983]  to  calculate 
the  effects  of  small  scale  ionization  on  a  propagating  wide  bandwidth  signal.  The  results  from 
averaging  many  Monte  Carlo  realizations  are  compared  to  the  Phase-Screen  Diffraction  Method 
(PDM)  of  Nickisch  [1992],  This  is  part  of  an  effort  to  validate  the  new  Matlab  implementation 
of  the  MPS  code. 

To  obtain  the  solution  for  wide  bandwidth  propagation,  the  MPS  code  is  applied  to  a 
spectrum  of  discrete  frequencies  over  the  bandwidth  of  the  signal  of  interest.  Note  that  the  MPS 
realizations  provide,  in  engineering  terminology,  the  transfer  function  of  the  ionospheric 
propagation  channel  at  a  specific  transmission  frequency. 

A  useful  time  domain  signal  to  use  for  illustrative  purposes  is  the  triangular  function  given 
by  p(t)  =  1  -  |x|/Tc  for  -Tc  <  x  <  Tc.  The  Fourier  transform  is  P(f)  =  Tc  sin2  (7ifTc)/(  7ilTc)2.  This 
function  has  nulls  at  frequencies  that  are  multiples  of  2/Tc.  To  represent  this  function 
adequately,  we  use  MPS  frequencies  ranging  over  -2n/Tc  <  f  <  2n/Tc  where  n  is  an  optional 
integer,  usually  1  or  2. 

The  phase  screen  diffraction  method  (PDM)  [Nickisch,  1992]  technique  uses  split-step 
Fourier  transforms  to  solve  the  equation  for  the  mutual  coherence  function  under  the  strong 
scattering  assumption.  This  assumption  is  not  a  limitation  of  the  MPS  calculation.  The  PDM 
calculation  is  three-dimensional,  allows  for  phase  screens  that  move  in  time  with  different 
velocities,  is  extremely  fast  with  a  minimum  of  computer  resources  required,  and  is  an  excellent 
choice  for  comparison  to  MPS  calculations  for  strong  scattering  cases. 

In  the  following,  we  apply  the  MPS  propagation  simulation  to  calculate  some  important 
features  of  signals  that  propagate  through  structured  regions  of  the  ionosphere.  First  consider  the 
example  of  propagation  through  a  single  cloud  of  structured  ionization.  A  similar  example  was 
considered  in  Knepp  [1983],  Figure  D14  shows  an  example  realization  of  the  single  phase 
screen  used  to  represent  the  ionosphere. 

The  phase  screen  representation  here  consists  of  a  combination  of  two  major  features;  a 
large  Gaussian  lens  (to  model  the  gross  cloud)  with  a  half-size  of  5  km,  a  peak  phase  of  200 
radians  (which  gives  a  gross  propagation  delay  of  td  =  200/27ifc  =  0.32  microsec)  and 
superimposed  structured  ionization.  The  large  cloud  is  centered  in  the  MPS  calculation  grid  that 
consists  of  524,288  points  and  whose  length  is  50  km.  The  structured  portion  of  the  phase  of  the 
cloud  is  represented  with  a  power-law  power  spectral  density  (PSD)  with  q 3  slope  (which 
corresponds  to  a  q'4  for  the  3D  Ne  spectrum),  an  outer  scale  of  390  meters,  an  inner  scale  of  10 
meters,  and  a  phase  standard  deviation  of  4  radians.  These  values  are  not  realistic  for  the  natural 
ionosphere,  but  are  chosen  for  the  purpose  of  this  example.  A  total  of  64  spectral  components 
(frequencies)  are  calculated  over  a  total  bandwidth  of  20  MHz  centered  at  a  carrier  frequency  of 
fc  =  100  MHz.  The  propagation  geometry  consists  of  this  single  phase  screen  to  represent  the 
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ionization.  The  plane  wave  signal,  p(x),  propagates  a  distance  of  100  km  to  the  receiver.  The 
value  of  Tc  is  about  100  nsec. 


Delay  ( ju  sec )  Delay  (ju  sec) 

Figure  D15.  Time-domain  signal  after  prop-  Figure  D16.  Close  up  of  the  central  portion  of 
agation  100  km  past  the  cloud  of  ionization,  the  time-domain  signal  from  15  to  35  km  along 
Distance  along  the  MPS  grid  increases  from  top  the  MPS  grid, 
to  bottom. 


Figure  D14  shows  the  phase  of  the  single  phase  screen  as  a  function  of  distance  along  the 
numerical  MPS  grid.  Figure  D15  shows  a  realization  of  the  resulting  time  domain  waveform  as 
a  function  of  distance  along  the  MPS  grid.  In  Figure  D15,  100  plots  of  amplitude  versus  time  are 
shown  selected  uniformly  from  the  entire  MPS  grid  of  2 19  points.  At  the  edges  of  the  grid,  at  the 
top  and  bottom  of  the  figure,  the  triangular  pulses  are  not  affected.  The  amplitude  is  near  unity 
and  the  pulses  are  undistorted. 
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In  the  middle,  the  delay  is  greatest  where  the  ionization  structure  is  highest.  At  the  edges  of 
the  smooth  Gaussian,  there  is  also  some  focusing  where  the  signal  levels  are  noticeably  larger 
than  unity.  The  smooth  Gaussian  ionization  acts  like  a  defocusing  lens  with  a  minimum  near  the 
center  and  peaks  near  the  two  edges  (around  20  and  30  km).  In  the  center  of  the  grid  there  is 
significant  time  delay  jitter  which  is  caused  by  random  angular  scattering  that  is  a  function  of  the 
ionization  and  the  transmission  frequency.  Figure  D16  shows  a  close  up  of  the  interior  portion 
of  the  MPS  grid  from  15  to  35  km. 

The  MPS  propagation  code  contributes  usefully  to  the  analysis  of  the  effects  of  the 
ionosphere  on  propagating  wideband  signals.  For  wideband  signals,  representation  of  the 
propagating  signal  consists  of  a  complex  valued  function  of  both  time  and  frequency.  For  a 
complete  signal  representation,  one  essentially  needs,  at  a  single  instant,  the  amplitude  and  phase 
of  the  signal  at  each  spectral  component  within  the  system  bandwidth.  Useful  tools  to  help 
understand  the  impact  of  the  ionosphere  on  wide  bandwidth  transmission  systems  are  the 
impulse  response  function  and  its  Fourier  transform,  the  ionospheric  transfer  function.  For  wide 
bandwidth  signals,  the  MPS  code  directly  calculates  realizations  of  the  complex  I&Q  signal 
components  across  the  signal  spectrum.  This  calculation  gives  the  ionospheric  transfer  function. 
This  frequency-domain  function  is  then  multiplied  by  the  transfer  function  of  the  receiver  filters; 
then,  the  product  is  Fourier  transformed  to  obtain  the  time  domain  signal  that  is  input  into  the 
receiver  of  a  radar  or  communications  system. 

Figure  D17  shows  an  example  of  the  transfer  function  and  impulse  response  from  this  MPS 
calculation.  Simultaneous  portrayal  of  the  magnitude  of  the  transfer  function  and  that  of  the 
impulse  response  function  is  very  useful  in  understanding  the  effects  of  the  ionosphere  on  the 
propagating  wide  band  signal.  One  is  able  to  associate  effects  in  the  time  and  frequency  domains 
and  to  obtain  an  understanding  of  the  spatial  variation  of  the  signal  fluctuations.  It  is  also 
possible  to  employ  the  MPS  code  to  obtain  other  useful  statistical  averages  of  the  propagation 
channel.  Here,  the  MPS  realizations  are  analyzed  to  calculate  the  mutual  coherence  function 
T(x,f).  This  is  accomplished  by  Fourier  techniques  involving  multiple  crosscorrelations  of  each 
frequency  component  in  a  multi-frequency  realization  with  the  realization  at  the  carrier 
frequency.  This  result  is  then  Fourier  transformed  to  obtain  the  generalized  power  spectrum 
S(q,r). 

Figure  D18  shows  S(q,x)  for  another  example  of  an  MPS  calculation.  For  this  calculation, 
the  transmission  frequency  is  300  MHz  and  the  MPS  code  propagates  128  frequency  components 
over  a  total  two-sided  bandwidth  of  40  MHz  centered  at  a  carrier  frequency  of  300  MHz.  The 
propagation  geometry  consists  of  a  thick  scattering  layer  420  km  in  thickness,  which  is  modeled 

1  7 

using  10  equally  spaced  phase  screens.  The  MPS  grid  has  a  length  of  100  km  and  consists  of  2 
points.  The  signal  then  propagates  an  additional  420  km  to  the  receiver  plane.  The  total  rms 
phase  of  the  thick  layer  is  170  radians.  A  q 3  PSD  is  used  to  represent  the  phase  fluctuation  with 
an  outer  scale  of  5  km  and  inner  scale  of  10  m.  To  obtain  the  averages,  the  code  is  exercised  50 
times  using  a  different  random  seed  each  time. 

Figure  D19  shows  the  phase  screen  diffraction  method  (PDM)  [Nickisch,  1992]  solution  for 
the  same  propagation  geometry  as  that  used  to  create  Figure  D18.  The  two  figures  are  quite 
similar  and  validate  the  MPS  calculation. 
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Figure  D17.  Example  of  the  transfer  function  and  impulse  response  function  at  a  point  along  the  MPS  grid. 
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Figure  D18.  MPS  measurement  of  the  generalized  power  Figure  D19.  PDM  calculation  of  the  generalized 
spectrum  obtained  as  the  average  of  50  realizations.  power  spectrum. 

We  have  presented  a  comparison  of  the  multiple  phase  screen  results  for  the  generalized 
power  spectrum  S(q,x)  with  results  from  the  phase  screen  diffraction  method  (PDM).  The 
agreement  is  very  good  and  we  have  validated  the  Matlab  implementation  of  the  MPS  code.  In 
addition,  we  have  demonstrated  a  calculation  that  may  be  of  specific  interest  to  space  based  radar 
where  both  dispersion  and  diffraction  is  important. 
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Appendix  E 

Applications  of  GPS  TEC  Analysis  Tools 

Tools  developed  by  NWRA  under  previous  contract  efforts  and  extended  under  the  current 
contract  were  used  to  address  various  topics  of  interest  of  the  USAF.  In  this  Appendix  we 
describe  these  efforts.  In  particular,  we  supplied  ray-tracing  expertise  and  code,  and  we  provided 
analyses  based  on  our  ionospheric  inversion  algorithms  (GPSII  and  CREDO,  described  below). 
The  applications  included  AFRL  testing  of  various  ionospheric  models,  support  for  the  MDA 
Early  Launch  Detection  and  Tracking  program,  demonstration  of  potential  performance 
enhancements  for  a  specific  system  of  interest,  analysis  for  a  proposed  Next  Generation  OTH 
radar  system,  and  ionospheric  calibration  testing  for  Upgraded  Early  Warning  Radar. 

El.  Ionospheric  modeling  and  Ray  Tracing 

NWRA  delivered  ray-tracing  software  and  several  post-processing  routines  to  2nd 
Lieutenant  Shayne  Aune  of  AFIT  at  Wright-Patterson  AFB.  Lt.  Aune  was  preparing  a  master’s 
thesis  in  which  he  was  to  demonstrate  the  effect  of  the  ambient  ionosphere  upon  radio  wave 
propagation.  We  gave  him  an  executable  version  of  the  Jones-Stephenson  ray-tracing  program 
that  we  have  been  using  and  refining  for  our  own  use  for  many  years.  We  gave  instructions  in 
the  use  of  this  code,  and  acted  as  “tech  support”  for  the  questions  that  came  up  in  using  it. 

Among  the  output  generated  by  the  ray-tracing  program  are  binary  files  containing  ray 
trajectories.  We  developed  and  provided  Matlab  scripts  that  could  display  these  ray  trajectories 
in  an  easily  understood  format.  Figure  El  illustrates  one  example:  a  fan  of  rays  has  been  shot  at 
several  different  initial  elevation  angles  and  a  single  azimuth.  We  see  that  rays  with  low 
elevation  angles  are  reflected  back  to  ground,  while  higher-angle  rays  penetrate  the  ionosphere. 
The  color  contours  show  the  electron  density  of  the  ambient  ionosphere  in  the  same  plane  as  that 
of  the  ray  fan. 
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Figure  El.  A  fan  of  ray  trajectories,  with  the  ambient  ionosphere  shown. 

We  also  provided  a  version  of  our  executable  program  MakelRItable.  This  program  acts  as 
a  front  end  to  the  IRI-2001  model  of  the  ionosphere.  The  user  inputs  his  choice  of  a  latitude  and 
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longitude  grid,  sunspot  number,  and  time  and  date,  and  the  program  creates  a  three-dimensional 
table  of  ionospheric  electron  densities.  The  table  is  in  a  format  that  our  ray-tracing  program  can 
read.  The  color  contours  in  the  above  figure  show  an  example  of  MakelRItable  output. 

Program  MakelRItable  was  delivered  to  AFRL/RVBXP.  The  intent  was  to  produce  tables 
of  ionospheric  density  that  could  be  used  for  ray-tracing  studies  to  compare  the  results  of  the 
predicted  IRI  model  ionospheres  to  other  models,  such  as  GAIM.  AFRL  has  also  been 
examining  the  IRI  predictions  for  other  purposes.  They  noticed  a  peculiarity  in  the  dependence 
upon  sunspot  number  which  revealed  a  bug  in  our  MakelRItable  front  end.  This  bug  was  fixed, 
and  the  improved  version  of  MakelRItable  was  delivered  to  AFRL.  Issues  regarding  having  our 
ray  tracing  codes  run  successfully  at  AFRL  were  resolved. 

NWRA  provided  AFRL  with  our  GPS  Ionospheric  Inversion  code,  GPSII  [Fridman,  et  al., 
2006;  Nickisch,  et  al.,  2006],  AFRL,  with  our  assistance,  applied  this  code  in  a  Direction 
Finding  (DF)  test.  The  GPSII  algorithm  was  refined  and  tested,  and  a  User’s  Manual  was 
prepared  to  allow  users  at  AFRL  to  run  GPSII.  We  consolidated  various  subroutines  that  had 
appeared  in  slightly  modified  forms  between  our  CREDO  (referenced  below)  and  GPSII  codes; 
we  made  single  versions  that  are  compatible  with  both  of  these  codes  to  maintain  version  control. 

We  extended  the  GPSII  algorithm  to  incorporate  relative  TEC  data  from  LEO  satellites. 
Any  number  of  LEO  TEC  receivers  may  now  be  processed  simultaneously  within  GPSII 
together  with  data  from  GPS  receivers  and  vertical  sounders.  The  GPSII  system  automatically 
estimates  data  noise  using  sliding  windows.  Presently  the  system  reads  LEO  satellite  data  from 
FORTRAN  binary  files.  We  created  Matlab  scripts  for  converting  LEO  satellite  data  from  the 
NetCDF  format  to  FORTRAN  binary  files.  We  tested  the  operation  of  the  algorithm  by 
simultaneously  processing  data  from  19  GPS  receivers  and  2  LEO  receivers.  GPSII  is  still 
capable  of  working  in  real  time  on  a  2GHz  Pentium  PC.  We  also  incorporated  a  more  general 
RINEX  data  reading  capability  and  the  option  to  use  Cl  data  when  PI  data  is  unavailable.  The 
GPSII  user's  manual  was  updated  to  reflect  the  addition  of  the  LEO  capability. 

We  delivered  the  new  version  of  GPSII  with  the  capability  of  incorporating  beacon  data 
from  LEO  satellites.  The  GPSII  User’s  Guide,  updated  for  the  LEO  capability,  was  delivered  as 
well.  Assistance  on  GPSII  operation  was  provided,  on  several  occasions  to  Lt.  Curtis  Baragona. 
Assistance  on  both  GPSII  and  RAYTRACE  operation  was  provided  to  Lt.  Joshua  Warner.  A 
significant  upgrade  to  GPSII  was  also  delivered.  In  this  upgrade,  certain  memory-intensive 
matrix  manipulations  were  performed  in  a  new  way  that  does  not  require  the  full  expansion  of 
sparse  matrices.  It  makes  much  more  efficient  use  of  allocated  memory  than  prior  versions, 
resulting  in  a  much  smaller  memory  footprint  and  a  reduced  chance  of  crashing  due  to  a 
fragmented  memory  heap.  It  is  also  significantly  faster  for  large  data  sets  that  would  have 
stressed  the  older  version.  For  example,  in  a  case  that  involved  26  RINEX  receivers  collecting 
data  over  4  days,  it  took  almost  22  hours  to  process  with  the  old  executable;  the  same  case 
finished  in  10  hours  with  the  new  executable. 

NWRA  also  supported  a  MIMO  experimental  campaign  conducted  jointly  by  MIT  Lincoln 
Laboratory  and  the  Australian  DSTO  by  providing  predictions  of  propagation  conditions  for  the 
planned  campaign  time  and  location.  MIMO  (Multiple  Input  Multiple  Output),  as  applied  to 
OTH  radar,  refers  to  the  technique  of  transmitting  orthogonal  waveforms  from  the  OTH 
transmitter  elements,  allowing  transmit  beams  to  be  formed  after  signal  reception.  Experimental 
campaign  planning  required  knowledge  of  HF  skywave  propagation  mode  elevation  separations 
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to  ensure  that  the  antenna  apertures  were  properly  configured  to  be  able  to  isolate  the  modes 
using  MIMO.  Using  the  FIRIC  ionosphere  model  together  with  ray- tracing  using  our  ray¬ 
homing  algorithm,  we  provided  the  climatologically-predicted  propagation  mode  details  as  a 
function  of  radar  operating  frequency.  An  example  is  shown  in  Figure  E2.. 


OBI  file:  obi3d  14  11  2007  0400.out 
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Figure  E2.  Predicted  propagation  mode  elevations  at  the  transmitter  for  the  planned  MIMO 
experimental  link  from  JORN  Laverton  to  Broome,  Australia. 

We  also  supported  AFRL  by  providing  examples  of  output  from  our  ionospheric  inversion 
code  CREDO  (referenced  below)  for  briefings.  Using  ROTHR  data  for  two  24-hour  periods 
corresponding  to  low  and  high  sunspot  number,  we  created  ionosphere  models  for  each  24-hour 
period  (with  12-min  spacing)  and  made  movies  of  critical  plasma  frequency  over  the 
geographical  area  of  ROTHR  coverage.  An  example  of  one  of  the  movie  frames  is  shown  in 
Figure  E3. 

E2.  Early  Launch  Detection  and  Tracking 

Early  Launch  Detection  and  Tracking  (ELDT)  is  a  Missile  Defense  Agency  program  with 
the  goal  of  providing  a  means  to  detect  ballistic  missiles  immediately  after  launch.  SkyLOS  is  a 
series  of  experiments  supporting  ELDT  wherein  forward-based  receivers  are  used  to  detect  HF 
reflections  from  missile  launches  that  are  illuminated  by  a  remote  transmitter. 
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Critical  Frequency  (MHz),  26  Sep  1997, 18:55  EST 


Figure  E3.  Contours  of  ionospheric  critical  frequency  obtained  by  CREDO  inversion  of  ROTHR 

Virginia  sounding  data. 

Our  approach  was  to  use  CREDO  (Coordinate  Registration  Enhancement  by  Dynamic 
Optimization),  a  software  package  for  the  generation  of  ionosphere  models  from  OTH  radar 
soundings  [Fridman  and  Nickisch,  2001],  With  CREDO,  we  read  in  the  VI  and  WSBI  of  each 
launch  time,  generated  corresponding  ionosphere  models,  and  synthesized  the  OIs  corresponding 
to  the  DSTO  measurements.  The  synthetic  OIs  were  overlain  on  the  measured  OIs  to  confirm 
the  adequacy  of  the  ionospheric  solutions.  A  Matlab  translation  script  was  generated  to  convert 
the  CREDO  output  model  to  the  data  format  required  by  VTRPE,  and  the  ionosphere  models 
were  delivered  to  Jim  McGraw  (Photon  Research  Associates;  PRA)  for  VTRPE  analysis.  One  of 
the  two  launch  cases  is  shown  in  Figure  E4,  for  14  NOV  2005  at  2030  UT. 

Figure  E4  displays  the  WARF-measured  VI.  The  synthetic  VI  generated  from  the  CREDO 
solution  ionosphere  model  is  overlain,  indicating  that  a  good  matching  ionosphere  model  was 
obtained.  Figure  E5  displays  the  WARF  Wide  Sweep  Backscatter  Ionogram  (WSBI).  Overlain 
are  the  leading  edges  of  major  propagation  modes  from  the  fitted  CREDO  ionosphere  model. 
Also  overlain  is  the  synthetic  OI  from  the  WARF  transmit  site  to  WSMR.  Note  that  the  endpoint 
of  the  WSBI  is  at  the  WARF  receive  site  whereas  the  OI  is  transmitted  from  the  WARF  transmit 
site,  which  is  why  the  nose  of  the  OI  does  not  exactly  align  with  the  WSBI  leading  edge. 
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Figure  E4.  WARF  VI  with  overlain  CREDO  solution. 


Figure  E5.  WARF  WSBI  with  overlain  CREDO  propagation  mode  leading  edges  and  synthetic  OI. 
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Figure  E6  shows  the  DSTO-measured  01  with  the  synthetic  CREDO-based  01  overlain. 
Note  the  obvious  Traveling  Ionospheric  Disturbance  (TID)  activity  evident  in  the  measurement, 
particularly  as  a  ripple  in  the  F2  high-ray  trace  and  again  at  the  01  F2  nose.  These  cannot  be 
modeled  in  the  CREDO  ionosphere  since  the  WSBI  sounding  bandwidth  and  antenna  aperture 
size  do  not  provide  sufficient  resolution.  Aside  from  this,  the  match  is  quite  good. 


Figure  E6.  DSTO-measured  OI  and  overlain  CREDO  synthetic  OI. 

Figure  E7  shows  the  CREDO  ionosphere  model,  plotted  as  plasma  frequency  in  MHz. 
Figure  E8  displays  the  signal  strength  (relative  to  free-space  propagation)  of  the  full  wave 
VTRPE  code  using  the  CREDO  ionosphere  model.  The  transmit  antenna  elevation  profile  was 
assumed  to  be  that  of  a  vertical  dipole.  Near  the  bottom  right  hand  comer  of  the  plot  is  the 
trajectory  of  the  missile  launch.  Ongoing  analysis  is  being  carried  out  at  PRA  to  determine  how 
well  the  CREDO  ionosphere  model  predicts  measured  values  of  SkyLOS  parameters  like  time- 
difference  of  arrival  and  frequency-difference  of  arrival. 

The  SkyLOS  IV  campaign  was  conducted  in  June  2006.  Again,  the  WARF  OTH  radar 
transmitter  in  California  was  used  to  illuminate  two  missile  launches  from  White  Sands  Missile 
Range  (WSMR).  For  the  previous  SkyLOS  III  campaign  discussed  above,  Vertical  Incidence 
(VI)  and  Wide-Sweep  Backscatter  Ionograms  (WSBIs)  were  collected  by  WARF  to  assist  in 
modeling  the  ionosphere.  These  were  unavailable  for  SkyLOS  IV.  However,  Oblique 
Ionograms  (OIs)  were  collected  by  the  Australian  DSTO  on  a  path  from  the  WARF  transmit  site 
to  WSMR.  We  again  were  asked  to  provide  ionospheric  models  based  on  these  OIs  that  could  be 
used  in  the  full-wave  propagation  code  VTRPE. 

As  with  SkyLOS  III,  we  used  CREDO  for  generating  the  ionosphere  models,  but  were 
hampered  by  the  fact  that  CREDO  takes  as  input  VI  and  WSBI  soundings,  but  not  OIs,  and  VI 
and  WSBIs  from  WARF  were  not  available  for  SkyLOS  IV  (the  WARF  transmit  site  had  been 
dismantled  by  this  time).  We  resorted  to  converting  the  OIs  to  equivalent  Vis  using  the  well- 
known  Equivalence  theorems  of  ray  tracing.  These  we  supplemented  with  Vis  at  Pt.  Arguello 
and  Dyess  AFB,  which  we  obtained  from  DIDBase,  an  archive  of  vertical  soundings  maintained 
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Figure  E7.  CREDO  ionosphere  model. 
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Figure  E8.  VTRPE  full-wave  signal  strength  plot  in  CREDO  ionosphere. 
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by  the  University  of  Lowell  Center  for  Atmospheric  Research.  We  generated  corresponding 
ionosphere  models,  and  synthesized  the  OIs  corresponding  to  the  DSTO  measurements.  The 
synthetic  OIs  were  overlaid  on  the  measured  OIs  to  confirm  the  adequacy  of  the  ionospheric 
solutions.  The  CREDO  output  model  was  converted  to  the  data  format  required  by  VTRPE,  and 
the  ionosphere  models  were  delivered  to  Jim  McGraw  (Photon  Research  Associates;  PRA)  for 
VTRPE  analysis.  We  also  developed  and  supplied  to  PRA,  a  Hamiltonian  magnetoionic  3D  ray¬ 
tracing  capability  for  their  use  in  the  SkyLOS  data  analysis.  We  assisted  them  in  the  proper  use 
of  the  software. 

NWRA  assisted  in  planning  for  the  SkyLOS  V  campaign.  In  SkyLOS  V,  the  Relocatable 
Over-the-Horizon  Radar  (ROTHR)  transmitter  in  Virginia  was  used  to  illuminate  a  launch  of  the 
space  shuttle  Endeavor  from  Cape  Canaveral,  FL.  A  fundamental  issue  with  planning  for 
SkyLOS  V  was  determining  what  the  area  coverage  would  be  from  the  Virginia  ROTHR  during 
the  months  May  -  September  2007.  Our  solution  was  to  generate  FIRIC  model  ionospheres  for 
several  dates  within  that  period,  using  an  assumed  value  for  the  sunspot  number.  Based  upon  the 
estimates  available  from  NOAA,  we  estimated  a  most-likely  sunspot  number  of  about  12  for  this 
time  period.  A  few  models  were  created  with  sunspot  numbers  6  and  25,  to  test  sensitivity  to  this 
parameter;  we  found  that  this  variation  had  little  effect  upon  the  expected  coverage  area  for  any 
given  frequency.  We  then  ran  our  code  “Raytrace”  to  follow  the  trajectories  of  rays  launched 
into  these  ionosphere  models.  By  launching  fans  of  rays  at  numerous  azimuths  between  the 
limits  of  ROTHR’ s  coverage,  we  were  able  to  estimate  coverage  at  various  frequencies  by  noting 
the  landing  locations  of  those  rays  that  were  reflected  back  to  earth.  We  also  developed  a 
method  to  estimate  the  areal  power  density  represented  by  these  rays  (i.e.,  received  power  per 
unit  area  at  ground  level)  and  display  them  upon  a  geographic  map.  An  example  is  shown  in 
Figure  E9. 


Power  density  at  8.0  MHz,  03  Aug  2007,  UT  00:00,  SSN  =  12  dBm 
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Figure  E9.  Sample  coverage  map  generated  for  SkyLOS  V. 
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Similar  maps  were  generated  for  24-hour  time  spans  on  several  dates  (August  3,  June  28, 
and  September  7).  Frequencies  were  chosen  to  ensure  that  ample  power  reached  the  launch  site 
at  Kennedy  Spaceflight  Center  (shown  as  a  blue  star  on  the  figure  above). 

The  possible  effects  and  intensity  of  Sporadic-E  (Es)  for  the  dates  of  interest  were  also 
estimated.  A  survey  of  literature  and  past  ROTHR  ionogram  data  allowed  us  to  make 
probabilistic  predictions  of  Es  intensity  and  frequency  of  occurrence;  it  appeared  that  Es  would 
likely  be  somewhat  less  of  a  problem  than  it  was  for  the  SkyLOS  IV  experiments.  The  effect  of 
sporadic-E  on  individual  measurements  is  of  course  unpredictable,  but  we  can  roughly  estimate 
the  effects  of  Es  upon  coverage  areas.  We  re-ran  Raytrace  with  a  very  narrow,  strong,  Es-like 
layer  replacing  the  standard  E  layer.  The  almost  mirror-like  Es  reflection  means  that  distances 
much  closer  to  the  antenna  site  can  be  reached  at  any  given  frequency.  We  then  combined  the 
power  density  for  this  simulation  with  the  results  of  an  Es-free  simulation  to  get  a  coverage  map 
roughly  representing  50%  opacity  of  the  layer.  An  example  of  the  resulting  coverage  map  is 
shown  in  Figure  E10. 


Figure  E10.  Coverage  map  combining  power  from  Es  and  normal  propagation. 

Using  CREDO,  we  also  developed  a  similar  set  of  ionospheric  HF  propagation  predictions 
for  the  SkyLOS  V  campaign  using  actual  ROTHR-VA  soundings  for  times  of  interest  during  the 
month  preceding  the  actual  launch  on  August  8,  2007.  These  were  used  by  Photon  Research 
Associates  as  the  propagation  environment  for  a  VTRPE  full  wave  analysis.  Following  the 
actual  launch  and  data  collect  by  the  SkyLOS  team,  we  received  the  ROTHR  soundings  and 
oblique  ionograms  collected  by  the  Australian  DSTO  during  the  launch.  These  were  used  by  us 
to  develop  models  of  the  real  ionosphere  in  effect  during  the  SkyLOS  Y  campaign,  and  were 
passed  to  PRA  for  VTRPE  analysis  and  subsequent  TDOA/FDOA  measured  trajectory  analysis. 
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Consideration  was  given  to  the  best  means  for  automatic  frequency  management  for 
SkyLOS.  Thus  far,  the  ionospheric  modeling  for  SkyLOS  has  been  performed  using  CREDO, 
which  uses  OTH  radar  soundings  to  determine  the  current  ionospheric  state,  and  which  can  be 
used  to  predict  optimal  frequencies  for  illumination  of  the  launch  site  by  a  remote  HF  transmitter 
(as  well  as  determining  the  HF  propagation  paths,  a  requirement  for  accurate  FDOA 
geolocation).  Ultimately,  however,  the  SkyLOS  concept  must  be  viable  for  situations  that  are 
devoid  of  OTH  radar  support.  We  have  suggested  that  the  ionospheric  modeling  capability  could 
be  performed  with  GPSII,  the  NWRA  GPS  Ionospheric  Inversion  capability.  GPSII  uses  GPS 
two-frequency  beacon  transmissions  to  model  the  current  ionospheric  state,  optionally 
augmented  with  vertical  sounding  data.  GPSII  is  a  good  candidate  for  SkyLOS  because  the 
SkyLOS  transmitter  and  all  receiver  sites  are  already  equipped  with  GPS  receivers  (for  timing 
purposes). 

E3.  CREDO  Application  to  a  Specific  System 

AFRL  provided  NWRA  with  sounding  data  from  a  specific  system-of-interest.  We 
developed  the  capability  to  read  and  display  these  data  in  Ionoview,  the  graphical  user  interface 
for  our  CREDO  capability  (Coordinate  Registration  Enhancement  by  Dynamic  Optimization). 
We  made  modifications  to  our  CREDO  code  in  order  to  be  able  to  develop  three-dimensional 
ionosphere  models  from  these  sounding  data.  These  were  used  to  perform  a  demonstration  of 
CREDO’S  capability  of  providing  high-fidelity  coordinate  registration  for  this  system.  All  of  the 
ionograms  that  we  received  were  examined.  A  time  sequence  of  ionosphere  models  was  made 
for  the  relevant  geographic  region  covering  the  three-hour  period  for  which  ionogram  data  were 
available. 

CREDO  was  extended  to  allow  the  generation  of  three-dimensional  ionosphere  models, 
even  when  no  Backscatter  Ionogram  (BI)  data  are  available,  or  when  only  single  sector  Bis  are 
available  (in  addition  to  vertical  soundings).  This  makes  CREDO  useful  in  regions  where 
vertical  sounding  data  are  available,  but  where  there  may  be  no  OTH  radar  to  provide  backscatter 
ionograms,  or  for  special  systems  with  only  a  single  BI  sector.  Previously,  CREDO  worked  with 
single-sector  WARE  radar  soundings,  but  provided  only  a  two-dimensional  (range-altitude) 
ionosphere  model  in  the  direction  of  the  single-sector  BI.  The  extended  version  of  CREDO  uses 
the  underlying  FIRIC  ionosphere  model  to  initiate  the  ionosphere  model  beyond  the  single 
sector,  expanding  to  nine  sectors  covering  90  degrees  of  azimuth.  That  initial  model  is  then 
adjusted  to  derive  an  ionosphere  model  that  matches  the  single-sector  BI  data  and  available 
vertical  sounding  data. 

We  applied  this  capability  to  data  from  this  system.  We  also  developed  new  code  for 
Coordinate  Registration  (CR)  of  targets  using  CREDO-generated  ionosphere  models.  The  CR 
routine  computes  fans  of  rays  from  the  transmitter  and  receiver  locations,  allowing  the 
construction  of  all  possible  two-way  propagation  modes  to  given  target  radar  slant  coordinates 
(delay  and  azimuth-of-arrival  at  the  receiver).  A  mode  association  algorithm  is  applied  that 
seeks  to  minimize  the  variance  of  inferred  target  ground  coordinates  as  weighted  by  computed 
mode  signal  strength.  The  procedure  was  applied  to  test  data  with  GPS  truth.  The  CREDO 
improvement  over  the  currently  implemented  method  was  documented  and  delivered  to  the 
customer. 
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E4.  Next  Generation  OTH  Radar  Analysis 

NORAD/NORTHCOM  is  promoting  the  development  of  a  next-generation  OTH  radar, 
dubbed  NexGen.  The  initial  phase  of  this  endeavor  is  to  develop  a  NexGen  Roadmap  that 
defines  the  requirements  for  such  a  system.  We  were  asked  to  support  the  definition  of  the 
Coordinate  Registration  function  of  the  NexGen  system. 

A  new  capability  of  the  NexGen  system  will  be  backscatter  soundings  received  on  a  two- 
dimensional  aperture,  allowing  some  measure  of  the  received  elevation  angles  of  the 
backscattered  propagation  modes.  We  developed  a  simulation  capability  for  generating  such 
backscatter  ionograms.  We  resurrected  the  Clutter  Effects  Model  (CLEM)  code  [Lauer,  et  al., 
1998],  which  models  OTH  radar  clutter,  including  Spread  Doppler  Clutter  (SDC)  caused  by 
naturally  occurring  small-scale  ionization  structure.  Ultimately,  we  hope  to  use  the  algorithms  in 
CLEM  as  the  basis  for  a  detailed  NexGen  simulator.  For  the  current  effort,  however,  we 
concentrated  on  the  simulation  of  Elevation  Backscatter  Ionograms  (EBIs).  Backscatter 
Ionograms  (Bis)  of  current  OTH  radars  can  provide  only  frequency-delay  information  about  the 
propagation  mode  environment.  Unfolding  this  information  in  receive  elevation  will  provide  a 
considerable  additional  constraint  on  inverse-processing  solutions  for  the  current  ionosphere 
model,  which  should  significantly  improve  Coordinate  Registration  accuracy.  Using  CLEM,  we 
were  able  to  synthesize  an  EBI  to  illustrate  this  extra  information.  Figure  Ell  displays  a 
synthetic  BI  created  using  CLEM,  but  displayed  in  the  traditional  frequency-delay  format.  In 
Figure  El 2,  the  view  is  rotated  to  expose  the  receive  elevation  axis.  Note  that  E-mode  returns 
are  separated  from  the  F-mode  returns  in  this  view  and  are  easily  distinguished. 


Figure  Ell.  Synthetic  Backscatter  Ionogram  in  traditional  frequency-delay  (slant  range)  format. 
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Figure  E12.  Synthetic  Elevation  Backscatter  Ionogram. 
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Appendix  F 

Example  Output  From  the  SCORE  Implementation 

This  appendix  provides  examples  of  graphical  displays  that  were  developed  under  this 
contract  effort  to  provide  tools  for  better  assessment  of  the  accuracy  of  GPS  observations  derived 
from  measurements  of  ionospheric  group-delay  and  phase-advance.  As  described  in  Section 
3.2.1.  of  this  report,  these  tools  were  developed  to  provide  better  assessment  of  the  performance 
of  the  SCORE  system  and,  once  it  is  ready,  for  the  SCORPION  system. 

Three  different  types  of  plots  are  generated  by  the  standard  GPS  TEC  processing: 

a.  Latitude-slice  plots  of  TEC  as  a  function  of  time  and  IPP  latitude  as  a  function  of  time. 
Typically,  more  than  one  latitude- slice  plot  is  generated  per  station-day.  Figure  FI  shows  an 
example  of  an  “all”  latitude  slice,  and  Figure  F2  is  an  example  of  a  “center”  latitude  slice 
showing  the  inclusion  of  the  average  TEC  and  +/-  one  standard  deviation  from  the  average 
(heavy  and  light  black  curves,  respectively). 

b.  Bias  comparison  plots,  plotting  the  biases  derived  in  the  SCORE  processing  against  the 
biases  from  the  AIUB  fdes.  Figure  F3  is  an  example  of  this  plot. 

c.  By-PRN  TEC  average  and  variation  plots.  This  is  a  plot,  by  satellite,  of  the  average 
difference  and  standard  deviation  between  the  TEC  for  that  satellite  and  the  mean  TEC 
calculated  from  the  “center”  latitude  slice.  Figure  F4  is  an  example  of  this  plot.  The  two 
horizontal  dashed  lines  are  the  average  standard  deviation  for  the  entire  data  set,  and  the 
vertical  bars  at  the  bottom  indicate  the  relative  number  of  points  by  PRN. 

We  were  able  to  use  TEC  observations  from  the  Jason- 1  satellite-based  altimeter  for 
comparisons  with  the  GPS  TEC  data  as  part  of  the  quality  control  and  uncertainty  assessment 
processes.  Two  different  types  of  plots  are  generated  in  the  GPS-Jason  comparison: 

a.  A  “composite”  plot  showing  a  variety  of  parameters  from  the  comparison.  Figure  F5 
shows  an  example  of  this  (in  geographic  coordinates).  The  red  triangles  are  the  sub-satellite 
latitude  and  longitude  for  the  Jason  points  used  in  the  comparison,  the  blue  circles  are  the 
(smoothed  and  resampled)  Jason- 1  TEC  values  as  a  function  of  latitude,  the  black  bars 
indicate  the  weighted  GPS  average  and  standard  deviation  TEC  as  a  function  of  latitude,  and 
the  green  diamonds  are  the  weighted  average  difference  between  the  Jason- 1  and  GPS  TEC. 
The  bias  file  used  to  create  the  GPS  TEC  data,  and  the  settings  for  the  parameters  used  to 
control  the  comparison  process  are  given  at  the  bottom-left  of  the  plot. 

b.  A  “multiple”  plot  showing  four  different  views  of  the  data  included  in  the  composite  plot. 
There  are  two  versions  of  this,  one  which  has  a  plot  of  the  azimuth  and  elevation  of  both  the 
GPS  and  Jason  TEC  points  and  one  which  has  a  plot  of  the  latitude  and  longitude  of  these 
points.  Figure  F6  is  an  example  of  the  first  type,  and  Figure  F7  is  an  example  of  the  second 
type.  These  plots  show  individual  data  points,  not  averages.  In  each  plot,  the  Jason  data  are 
shown  as  black  triangles  and  the  GPS  data  are  shown  as  colored  circles,  with  the  color 
indicating  which  satellite  the  data  came  from.  The  other  three  plots  on  each  figure  are,  from 
bottom  left  clockwise,  local  time  as  a  function  of  latitude,  TEC  as  a  function  of  latitude,  and 
TEC  as  a  function  of  local  time. 
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We  also  developed  a  set  of  after-the-fact  plots  providing  various  views  of  the  behavior  of  the 
satellite  and  receiver  biases  as  derived  from  both  the  raw  GPS  data  during  the  SCORE  analysis 
and  by  comparison  with  biases  obtained  from  the  AIUB. 

a.  A  pair  of  plots  for  a  specific  station  showing  the  variation  with  time  of  (1)  the  nominal 
receiver  bias  derived  from  fitting  the  SCORE-produced  biases  with  AIUB  biases,  and  (2)  the 
slope  of  the  linear  fit  derived  from  comparing  the  SCORE-produced  biases  with  AIUB 
biases  (the  slope  of  the  linear  fit  shown  in  Figure  F3).  Figure  F8  is  an  example  of  this  plot 
for  the  station  at  Neah  Bay,  WA.  Note  that  the  receiver  bias  plotted  is  the  difference 
between  each  individual  day’s  estimated  receiver  bias  and  the  mean  bias  for  the  full  set 
shown  in  the  plot.  Sudden  changes  in  the  character  of  these  plots,  such  as  the  jump  in  the 
receiver  bias  on  day  94,  can  be  used  to  watch  for  the  need  to  use  a  new  set  of  biases. 

b.  A  plot  for  a  specific  station  of  the  combined  satellite+receiver  bias  for  each  PRN  as  a 
function  of  time.  Figure  F9  is  an  example  of  this  plot  for  Neah  Bay,  WA,  for  the  same  time 
period  covered  by  the  data  in  Figure  F8.  As  in  the  receiver-bias  plot  in  Figure  F8,  shown 
plotted  is  the  difference  between  the  combined  bias  for  each  day  and  the  mean  bias  for  the 
entire  period  for  that  PRN.  A  sudden  change  in  all  PRNs,  such  as  occurs  at  day  94,  can 
indicate  a  change  in  the  receiver  bias,  while  a  change  in  only  one  PRN  may  indicate  either  a 
change  in  that  satellite’s  bias  or  a  change  in  the  local  propagation  environment  which  results 
in  a  change  in  the  multipath  characteristics  for  that  PRN. 
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VTEC  (TECU) 


Station:  stjo(47.596, 307.322) 

Latitude  band:  all  (All  Latitudes)  Date(UT):  20031217 

Elevation  Limit:  10.0 
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Figure  FI.  Sample  SCORE  “all”  latitude  slice  plot. 
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VTEC  (TECU) 


Station:  stjo(47.596, 307.322) 

Latitude  band:  center  (46.596  to  48.596)  Date(UT):  20031217 

Elevation  Limit:  10.0 


IPP  Solar  Local  Time 

Figure  F2.  Sample  SCORE  “center”  latitude  slice  plot. 
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Compare  SCORE  to  AIUB  Biases  (stjo  3351) 


Figure  F3.  Sample  SCORE  AIUB  bias  comparison  plot. 
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Station:  stjo(47.596, 307.322) 

Latitude  band:  (46.60,48.60)  Date(UT):  2003 1217 

Elevation  Limit:  10.00 


(Vertical  bars  indicate  relative  #obs/PRN) 
Figure  F4.  Sample  SCORE  by-PRN  quality-control  plot. 
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GPS/Jason-1  TEC  Comparison 


St.  John’s,  Newfoundland 
Lat:  47.596  MLat:  54.386 
Lon:  307.322  MLon:  30.673 
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File  Plotted:  teccomp_Rljo_20031218_1827.dat 
Bias  File  Used:  bias3352-stjo-std.out 

Comparison  Settings:  Lat  (3.0,l<0)deg  /  LT  (0,3 ,0,1)  hrs  /  UT  (3,0,3-0)  hrs 


Figure  F5.  Example  of  the  GPS-Jason  comparison  composite  plot. 
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Figure  F6.  Example  of  the  GPS-Jason  multiple  plot  (with  elevation-by  azimuth  option). 
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Figure  F7.  Example  of  the  GPS-Jason  multiple  plot  (with  latitude-by-longitude  option). 
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350km  IPP  Latitude  (deg) 


Slope  Receiver  Bias  -  Mean 


Station:  Neah  Bay,  WA 

(Geog:  48.30,235.38  I  Magn:  54.08,293.93) 


Figure  F8.  Example  of  an  individual  station  plot  of  receiver  bias  and  slope  as  derived  from  fitting  the 
station’s  biases  to  the  AIUB  biases. 
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Appendix  G 

Report  on  Research  by  Dr.  Spencer  Kuo 

The  following  report  was  provided  by  Dr.  Spencer  Kuo,  describing  work  he  completed  as  a 
consultant  to  NWRA  on  this  contract.  Note  that  further  details  of  the  reported  work  can  be  found 
in  the  publications  with  Dr.  Kuo  as  author  or  co-author  in  Section  4. 

Gl.  Introduction 

Several  prominent  issues  relevant  to  the  applications  of  the  HAARP  heating  facility  at 
Gakona,  Alaska  were  addressed.  Specifically,  the  topics  of  research  included  1)  Extremely/very 
low  frequency  (ELF/VLF)  wave  generation,  generation  mechanisms,  and  locating  the  electrojet 
source  region,  2)  precipitation  of  magnetospheric  energetic  electrons  by  whistler  waves,  which 
may  be  generated  by  the  HAARP  HF  heater,  3)  parametric  instabilities  and  cascades  of 
Langmuir  waves  and  upper  hybrid  waves,  which  produce  energetic  electrons  responsible  for 
airglow  enhancements,  and  4)  filamentation  instability  exciting  large-scale  field-aligned  density 
irregularities  and  geomagnetic  micropulsations.  Most  of  the  finished  works  have  been  published 
in  journal  articles. 

G2.  Extremely /Very  Low  Frequency  (ELF/VLF)  Wave  Generation 

The  conductivity  of  the  polar  electrojet  can  be  modulated  temporally,  through  Ohmic 
heating  by  power-modulated  HF  heating  waves,  to  set  up  the  ionospheric  antenna  current  for 
radiation.  The  dependency  of  the  radiation  intensity  of  this  virtual  ionospheric  antenna  on  the  HF 
wave  modulation  frequency  has  been  measured  in  HAARP  heating  experiments.  The  field 
intensity  of  the  ELF  wave  is  shown  to  be  affected  strongly  by  the  near  field  effect  (for  the 
modulation  frequency  fi  <  50  Hz)  and  cavity  resonance  effect  (for  fi  ~  2  kHz).  When  the  two 
effects  are  included  in  the  numerical  model,  the  experimentally  measured  frequency  dependency  of 
the  radiation  intensity,  shown  in  Figure  Gl,  can  be  reproduced  numerically.  On  the  other  hand,  the 
experimental  results  have  validated  the  numerical  model  on  ELF  wave  generation  by  the  HF 
modulated  electrojet. 


-HI 


-20 


hi  1  nr  io  J  to 1 

Frequency  (Ei/) 

Figure  Gl.  Frequency  dependency  of  ELF/VLF  radiation  intensity;  rectangular  wave  modulation  scheme  was 
used  in  the  experiments. 
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Experiments  using  the  intensity-modulated  HF  heating  waves  to  interact  with  electrojet 
currents  for  the  generation  of  VLF  waves  were  also  conducted.  An  unexpected  large  increasing 
rate  from  4  to  8  kHz  in  the  frequency  dependency  of  the  VLF  radiation  intensity,  shown  in 
Figure  Gl,  was  observed.  The  peak  value  at  8  kHz  was  intense  (about  7.5  dB  above  that  of  the  2 
kHz  signal  used  as  a  marker)  and  the  wave  intensity  from  5  to  17  kHz  appeared  to  be  abnormally 
high  (i.e.,  stronger  than  that  at  2  kHz).  In  the  experiments,  we  also  observed  the  enhancement  of 
spread-E  irregularities  at  electrojet  current  altitudes  due  to  the  amplitude-modulated  heater  wave. 
These  results  and  theoretical  analyses  suggest  that  temporally  modulated  electrojet  currents  mix 
with  heater  wave-excited  density  irregularities  to  form  whistler  mode  currents,  which  generate 
VLF  whistler  waves  directly  with  much  larger  intensities  and  better  directivity  than  a  Hertzian 
dipole  can. 


Figure  G2.  Radiation  intensity  ratios  (125/1205, 125/1325,  and  125/1625)  versus  altitude  of  ELF  source  in  the 
electrojet  current,  The  experimental  measurements  at  the  three  ratios  are  marked  by  on  the 
corresponding  plots. 

The  numerical  model  on  ELF  wave  generation  can  be  used  to  determine  the  altitude 
dependency  of  a  ratio  of  the  ELF  radiation  intensities  at  two  different  frequencies.  With  the  aid 
of  the  numerical  plot  on  the  altitude  dependency  of  an  intensity  ratio,  the  experimentally 
measured  single  value  of  a  similar  ratio  can  provide  unambiguous  information  on  the  height  of 
the  electrojet  contributing  to  the  ELF  radiation.  The  numerical  plots  of  three  independent 
intensity  ratios  (25  Hz/205  Hz,  25  Hz/325  Hz,  and  25  Hz/625  Hz)  are  presented  in  Figure  G2; 
these  plots  are  used  to  demonstrate  this  new  diagnostic  technique.  Each  ratio  of  the  measured 
ELF  radiation  intensities  is  also  marked  on  the  corresponding  plot,  with  the  same  pair  of 
frequencies,  by  a  symbol  The  location  of  this  mark  is  determined  by  the  intersection  point 
of  a  horizontal  dashed  line,  determined  by  the  measured  ratio,  and  the  plot.  The  altitude  of  the 
source  region  is  then  found  from  the  position  of  this  marked  point  on  the  horizontal  axis.  Ideally, 
the  three  plots  should  result  to  a  same  height.  In  reality,  various  errors  in  the  experimental 
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measurements  and  in  the  parameters  used  in  the  simulations  can  cause  some  variations.  Indeed, 
as  seen  in  Figure  G2,  the  determined  source  region  is  at  an  altitude  approximately  between  94 
Km  and  97  Km  from  the  ground  level.  However,  the  range  of  uncertainty  is  relatively  small, 
demonstrating  the  capacity  of  this  diagnostic  technique. 


G3.  Mitigation  of  Energetic  Electrons  in  the  Magnetosphere 

An  optimal  approach  reducing  the  population  of  MeV  electrons  in  the  magnetosphere  has 
been  explored.  It  invokes  a  double  resonance  condition;  whistler  wave  is  simultaneously  in 
cyclotron  resonance  with  keV  and  MeV  electrons.  Thus,  the  injected  whistler  waves  can  first  be 
amplified  by  the  background  keV  electrons  via  loss-cone  negative  mass  instability  to  become 
effective  in  precipitating  MeV  electrons  via  cyclotron  resonance  elevated  chaotic  scattering. 

The  theory  describing  the  amplification  of  whistler  wave  by  (tens  of  keV)  energetic 
electrons  in  the  magnetosphere  through  loss-cone  negative  mass  instability  is  formulated.  The 
numerical  results,  presented  in  Figure  G3a,  show  that  a  small  amplitude  whistler  wave  can  be 
amplified  by  more  than  25  dB.  The  amplification  factor  reduces  only  about  10  dB  with  a  30-dB 
increase  of  the  initial  wave  intensity.  As  the  wave  field  exceeds  the  threshold  for  the 
commencement  of  chaos,  cyclotron  resonance-enhanced  chaotic  scattering  can  precipitate  deeply 
trapped  MeV  electrons  into  loss  cones.  The  numerical  results,  presented  in  Figure  G3b, 
demonstrate  that  a  1.5  MeV  electron  can  be  scattered  from  an  initial  pitch  angle  of  86.5°  to  a 
pitch  angle  <50°  by  a  whistler  wave  with  the  magnetic  field  amplitude  of  0.08%  of  the 
background  magnetic  field,  which  is  about  20  times  smaller  than  that  without  invoking  cyclotron 
resonance.  This  percentage  converts  to  about  3100  pT  at  L=2,  and  200  pT  at  L=5,  which  reduce 
to  550  pT  and  36  pT,  respectively,  after  taking  advantage  of  the  15-dB  gain. 


Figure  G3.  (a)  Wave  amplitude  versus  time  for  four  different  incident  wave  intensities  and  (b)  temporal 
evolution  of  the  pitch  angle  of  a  1.5  MeV  electron  interacting  with  a  whistler  wave  at  Doppler  shifted 
cyclotron  resonance.  Plots  A  and  B  correspond  to  wave  magnetic  field  B1  =  7 ' 10-4  B0  and  8.1  '  10-4  B0 
cases,  respectively. 
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G4.  Parametric  Instabilities  and  Cascades  of  Langmuir  Waves 

Parametric  instability,  which  produces  cascade  spectra  of  HF  enhanced  plasma  lines 
(HFPLs)  in  HAARP  HF  heating  experiments,  is  studied.  In  each  cascade,  a  Langmuir  mother 
wave  decays  into  a  Langmuir  daughter  wave  and  an  ion  acoustic  wave;  The  Langmuir  daughter 
appears  to  be  an  HFPL  at  consecutively  decreasing  altitudes.  Thus,  the  cascade  lines  in  the 
HFPLs  spread  over  a  range  of  altitude,  and  the  threshold  pump  field  in  each  cascade  is  increased 
by  the  additional  propagation  loss  of  the  Langmuir  mother  wave.  The  dominant  factors,  which 
determine  the  number  of  cascade  lines  in  the  HFPLs,  include  the  ion-to-electron  temperature 
ratio,  T;/Te,  the  background  plasma  inhomogeneity  scale  length,  and  the  heating  wave  field 
intensity. 

The  ratio  Tc/T ,  (=  p)  depends  on  the  HF  heating  power  (i.e.,  oc  Ep ),  as  well  as  the  operating 
mode  of  the  HF  heater.  Assuming  that  the  HF  heater  is  operated  at  a  cw  mode,  electron  heating  is 
expected  and  is  modeled  by  a  function  r\(Ep)  =  1  +  1.5x[(fsp  -0.5)/4.5]2.  Using  this  heating 
model,  this  function  N(fsp;  p(fsp))  is  plotted  in  Figure  G4.  As  shown,  the  first  cascade  can  happen 
only  when  Ep  >  0.8  Y/m.  To  achieve  two  cascades,  Ep  has  to  exceed  1.8  Y/m.  As  the  heating 
power  increases,  the  electron  heating  gives  a  positive  feedback  to  further  increase  the  number  of 
cascade  lines  in  the  HFPLs.  At  the  maximum  heating  power  of  1  GW  ERP  (i.e.,  Ep  ~  3  V/m),  4 
discrete  cascade  lines  are  expected  to  appear  in  the  spectrum  of  HFPLs,  which  have  a  bandwidth 
exceeding  40  kHz.  These  lines  spread  in  altitude  over  a  0.5-km  range.  If  the  HF  heater  is 
operated  at  low  duty  cycle  pulsed  mode,  electron  heating  is  insignificant  even  at  a  very  high 
power  case.  Thus,  p  ~  1,  and  N0  <  3,  even  in  the  maximum  power  case  (i.e.,  Ep  ~  3  V/m). 


Figure  G4.  The  number  of  the  cascade  lines  NO  =  Integer  (N)  in  the  HFPLs  as  a  function  of  the  heating  wave 
amplitude  Ep,  for  h(Ep)  =  1  +  1.5 '[(Ep  -0.5)/4.5]2. 
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G5.  Filamentation  Instability  Exciting  Large-scale  Field-aligned  Density  Irregularites  and 
Geomagnetic  Micropulsations 

The  generation  of  large-scale  density  irregularities  by  the  filamentation  instability  near  the 
HF  reflection  height  was  evidenced  by  the  variation  of  the  reflection  height  of  the  sounding 
waves  in  the  frequency  region  around  the  HF  heater  frequency.  Significant  anomalous  absorption 
(i.e.,  reducing  the  backscattering)  of  the  sounding  waves  by  the  generated  density  irregularities 
caused  the  return  trace  of  the  sounding  waves  in  the  ionogram  to  disappear.  From  UTC  22:00  to 
23:00,  heater  wave  frequency  was  4.5  MHz.  As  shown  in  Figure  G5a,  the  return  of  the  sounding 
waves  recorded  in  the  ionogram  from  4.4  to  4.6  MHz  was  reduced  considerably.  After  UTC 
23:00,  the  heater  frequency  was  switched  to  2.7  MHz,  the  return  of  the  sounding  waves 
reappeared  in  the  ionogram.  Similarly,  as  shown  in  Figure  G5b,  the  return  trace  from  2.5  to  3 
MHz  in  the  ionogram  disappeared  during  the  period  from  23:00  to  01:00  UTC  when  the  heater 
was  run  at  2.7  MHz.  Again,  one  could  see  that  the  return  trace  reappeared  after  the  heater  was 
turned  off. 

Moreover,  as  shown  in  Figure  G5c,  the  magnetic  field  variations  in  D  and  Z  directions  were 
recorded  by  the  fluxgate  magnetometer  during  the  HF  heater  on  period.  These  results  presented 
in  Figure  G5,  showing  the  simultaneous  generation  of  large-scale  density  irregularities  and 
geomagnetic  field  fluctuations,  evidence  the  excitation  of  filamentation  instability. 


Figure  G5.  Composite  ionograms  showing  the  variation  of  the  reflection  height  of  the  sounding  waves  from 
(a)  4.4  to  4.6  MHz  and  (b)  2.5  to  3  MHz  every  15  minutes.  The  sounding  waves  appeared  to  be 
anomalously  absorbed  when  the  x  mode  HF  heater  at  4.5/2.7  MHz  was  turned  on;  (c)  magnetic  field 
variations  in  D  and  Z  directions  recorded  by  the  fluxgate  magnetometer. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


3D 

AEU 

AFB 

AFCRL 

AFIT 

AFRL 

AFWA 

AFS 

AK 

ALTAIR 

AMISR 

APTI 

ARPA 

AWS 

BI 

BMEWS 

cal/val 


Three  dimensional 
Antenna  Element  Unit 
Air  Force  Base 

Air  Force  Cambridge  Research  Laboratory  (now  AFRL) 

Air  Force  Institute  of  Technology 

Air  Force  Research  Laboratory 

Air  Force  Weather  Agency 

Air  Force  Station 

Alaska 

ARPA  Long-Range  Tracking  and  Instrumentation  Radar 
Advanced  Modular  Incoherent  Scatter  Radar 
Advanced  Power  Technology,  Inc. 

Advanced  Research  Projects  Agency 
Air  Weather  Service  (now  AFWA) 

Backscatter  Ionogram 

Ballistic  Missile  Early  Warning  Radar 

Calibration  and  Validation 


CLEM 

COSMOS 

CREDO 

DF 

DIDBase 

DMSP 

DSTO 

EBI 

ELDT 

ELF 

EOFs 

Es 

eV 


Clutter  Effects  Model 

Designator  for  a  certain  class  of  Russian  satellites 

Coordinate  Registration  Enhancement  by  Dynamic  Optimization 

Direction  Finding 

Digital  Ionogram  Data  Base 

Defense  Meteorological  Satellite  Program 

Defence  Science  Technology  Organisation 

Elevation  Backscatter  Ionogram 

Early  Launch  Detection  and  Tracking 

Extremely  Low  Frequency 

Empirical  Orthonormal  Functions 

Sporadic-E  ionosphere  layer 

Electron  volt 


f0F2  Critical  frequency  (O-mode)  of  the  F2  layer  (MHz) 

FAIM  Fully  Analytic  Ionosphere  Model 

FDOA  Frequency  Difference  of  Arrival 
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FIRIC 

FFT 

GAIM 

GAIM 

GIOS 

GPS 

GPSII 

HAARP 

HF 

I&Q 

IMS 

IPP 

IPT 

IRI 

IRI 

ITS 

JPL 

K 

kHz 

kW 

L-band 

LEO 

MDA 

MHz 

MIMO 

MPS 

NetCDF 

NexGen 

NIMS 

NNSS 

NOAA 

NORAD 

NORTHCOM 

NRE 


FAIM-IRI-Chapman  ionosphere  model 
Fast  Fourier  transform 

Global  Assimilative  Ionospheric  Model  (JPL  version) 

Global  Assimilation  of  Ionospheric  Measurements  (Utah  State  version) 

NWRA  GPS  Ionospheric  Observing  System  software 

Global  Positioning  System 

GPS  Ionospheric  Inversion 

High-frequency  Active  Auroral  Research  Program 

High  Frequency 

In-phase  and  quadrature-phase 

Ionospheric  Measuring  System 

Ionospheric  Penetration  Point 

Integrated  Product  Team 

International  Reference  Ionosphere  model 

Ionospheric  Research  Instrument  (located  at  HAARP  facility) 

Ionospheric  Tomography  System 

Jet  Propulsion  Laboratory 

Kelvin  (degrees) 

Kilohertz  (103  cycles/seconds) 

Kilowatt  (103  watts) 

Radio  frequency  band  covering  1.0  GHz  to  2.0  GHz  (nominal) 

Low  earth  orbit 
Missile  Defense  Agency 
Megahertz  (106  cycles/seconds) 

Multiple  Input  Multiple  Output 

Multiple  phase  screen 

Network  Common  Data  Format 

Next  Generation  OTH  radar 

Navy  Ionospheric  Monitoring  System 

Navy  Navigational  Satellite  System 

National  Oceanic  and  Atmospheric  Administration 

North  American  Aerospace  Defense  Command 

US  Northern  Command 

Northwest  Research  Engineering 
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NRL 

Naval  Research  Laboratory 

NTIA 

National  Telecommunications  and  Information  Administration 

NWRA 

Northwest  Research  Associates 

O&M 

Operations  and  Maintenance 

01 

Oblique  Ionogram 

ONR 

Office  of  Naval  Research 

OTH,  OTHR 

Over-the-Horizon  Radar 

PC 

Personal  Computer 

PCA 

Polar  Cap  Absorption 

PDM 

Phase  screen  diffraction  method 

PDR 

Powerful  Diagnostic  Radar 

PSD 

Power  spectral  density 

PEAK 

Propagation  Effects  Assessment  -  Kwajalein 

PMFIC 

Performance  Monitoring,  Fault  Isolation,  and  Calibration 

pot 

Potentiometer 

PRA 

Photon  Research  Associates 

PRN 

Pseudo-Random  Noise  (GPS  identification  signature) 

RCS 

Radar  cross  section 

RFP 

Request  for  Proposal 

RINEX 

Receiver  Independent  Exchange  (data  format) 

rms,  RMS 

root  mean  square 

ROTHR 

Relocatable  Over-the-Horizon  Radar 

ROTHR-VA 

Relocatable  Over-the-Horizon  Radar  located  in  Virginia 

S-band 

Radio  frequency  band  covering  2.0  GHz  to  4.0  GHz  (nominal) 

SAO 

Standard  Archiving  Output 

SCORE 

Self-Calibration  of  Range  Errors 

SCORPION 

SCORE  for  Plasmasphere  and  Ionosphere 

SDC 

Spread-Doppler  Clutter 

SEE 

Stimulated  Electromagnetic  Emission 

SkyLOS 

Skywave  Line  of  Sight 

SNR 

Signal-to-noise  ratio 

SOW 

Statement  of  Work 

SPAWAR 

Space  and  Naval  Warfare  Systems  Center 

sps 

samples  per  second 
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SRII 

SSIES 

SSUSI 

SSN 

SSPA 

TAC 

TDOA 

TEC 

TECU 

TLE 

UEWR 

UHF 

UK 

UML 

UPS 

UT 

VI 

VHF 

VLF 

VTEC 

VTRPE 

WARF 

WBMOD 

WDLS 

WSBI 

WSMR 

(Tj, 


SRI  International 

Special  Sensor  for  Ions,  Electrons,  and  Ions  (DMSP  sensor  package) 

Special  Sensor  Ultraviolet  Spectrographic  Imager  (DMSP  sensor  package) 

Sunspot  Number 

Solid  state  power  amplifier 

Technical  Advisory  Committee 

Time  Difference  of  Arrival 

Total  Electron  Content  (el/m2) 

TEC  unit  (l.OxlO16  el/m2) 

Two-Line  Element  satellite  ephemeris  set 

Upgraded  Early  Warning  Radar 

Ultra  High  Frequency  radio  band  (300  MHz  -  3  GHz) 

United  Kingdom 

University  of  Massachusetts  at  Lowell 
Uninterruptible  Power  Source 
Universal  Time 
Vertical  Incidence 

Very  High  Frequency  radio  band  (30  MHz  -  300  MHz) 

Very  Low  Frequency  radio  band  (3  kHz  -  30  kHz) 

Vertical  (or  equivalent  vertical)  TEC  (el/m2) 

Variable  Terrain  Radio  Parabolic  Equation 
Wide  Aperture  Research  Facility 

WideBand  MODel  (NWRA  ionospheric  scintillation  model) 

Weighted,  Damped,  Least-Squares 
Wide-Sweep  Backscatter  Ionogram 
White  Sands  Missile  Range 
Ionospheric  phase  scintillation  index  (radians) 
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